2009-12-01 28 views
36

malla de Numpy es muy útil para convertir dos vectores a una cuadrícula de coordenadas. ¿Cuál es la forma más fácil de extender esto a tres dimensiones? Por lo tanto, dados tres vectores x, y y z, construye matrices 3x3D (en lugar de matrices 2x2D) que se pueden usar como coordenadas.malla de malla Numpy en 3D

Respuesta

26

Aquí está el código fuente de meshgrid:

def meshgrid(x,y): 
    """ 
    Return coordinate matrices from two coordinate vectors. 

    Parameters 
    ---------- 
    x, y : ndarray 
     Two 1-D arrays representing the x and y coordinates of a grid. 

    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`. 

    See Also 
    -------- 
    index_tricks.mgrid : Construct a multi-dimensional "meshgrid" 
         using indexing notation. 
    index_tricks.ogrid : Construct an open multi-dimensional "meshgrid" 
         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]]) 

    `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) 

    """ 
    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) 

    y = y.reshape(numRows,1) 
    Y = y.repeat(numCols, axis=1) 
    return X, Y 

Es bastante sencillo de entender. Extendí el patrón a un número arbitrario de dimensiones, pero este código no está de ninguna manera optimizado (y tampoco está completamente verificado), pero obtienes lo que pagas. Creo que sirve:

def meshgrid2(*arrs): 
    arrs = tuple(reversed(arrs)) #edit 
    lens = map(len, arrs) 
    dim = len(arrs) 

    sz = 1 
    for s in lens: 
     sz*=s 

    ans = []  
    for i, arr in enumerate(arrs): 
     slc = [1]*dim 
     slc[i] = lens[i] 
     arr2 = asarray(arr).reshape(slc) 
     for j, sz in enumerate(lens): 
      if j!=i: 
       arr2 = arr2.repeat(sz, axis=j) 
     ans.append(arr2) 

    return tuple(ans) 
+1

En el caso de una malla de malla 3d, usando una muestra como la proporcionada en numpy doc para meshgrib, esto devolvería Z, Y, X en lugar de X, Y, Z. Reemplazando la declaración de devolución por 'return tuple (ans [:: - 1])' puede solucionar esto. – levesque

+0

@Paul si la longitud de la matriz x o y es larga, entonces el comando x.repeat() falla y envía un error de memoria. ¿Hay alguna forma de evitar este error? – Dalek

+0

@Dalek ¿Cuánto duran las matrices? ¿Podría ser que te quedas sin memoria? Por ejemplo, si hay 3 matrices, 4096 entradas cada una, y cada entrada contiene un doble (es decir, 8 bytes), entonces solo para las entradas que necesitamos (8 * 4 * 2 ** 10) ** 3 bytes = 2 ** 45 bytes = 32 * 2 ** 40 bytes = 32 TB de memoria, que obviamente es enorme. Espero no haberme equivocado aquí. –

5

Creo que lo que quiere es

X, Y, Z = numpy.mgrid[-10:10:100j, -10:10:100j, -10:10:100j] 

por ejemplo.

+0

Gracias, pero esto no es exactamente lo Necesito - meshgrid realmente usa los valores de los vectores para generar el arreglo 2D, y los valores pueden estar espaciados irregularmente. – astrofrog

7

¿Puede mostrarnos cómo está utilizando np.meshgrid? Existe una gran posibilidad de que realmente no necesite meshgrid porque la difusión numpy puede hacer lo mismo sin generar una matriz repetitiva.

Por ejemplo,

import numpy as np 

x=np.arange(2) 
y=np.arange(3) 
[X,Y] = np.meshgrid(x,y) 
S=X+Y 

print(S.shape) 
# (3, 2) 
# Note that meshgrid associates y with the 0-axis, and x with the 1-axis. 

print(S) 
# [[0 1] 
# [1 2] 
# [2 3]] 

s=np.empty((3,2)) 
print(s.shape) 
# (3, 2) 

# x.shape is (2,). 
# y.shape is (3,). 
# x's shape is broadcasted to (3,2) 
# y varies along the 0-axis, so to get its shape broadcasted, we first upgrade it to 
# have shape (3,1), using np.newaxis. Arrays of shape (3,1) can be broadcasted to 
# arrays of shape (3,2). 
s=x+y[:,np.newaxis] 
print(s) 
# [[0 1] 
# [1 2] 
# [2 3]] 

El punto es que S=X+Y puede y debe sustituirse por s=x+y[:,np.newaxis] porque este último no requiere (posiblemente grandes) arrays repetitivas a formar. También se generaliza a mayores dimensiones (más ejes) fácilmente. Simplemente agregue np.newaxis donde sea necesario para efectuar la transmisión según sea necesario.

Consulte http://www.scipy.org/EricsBroadcastingDoc para más información sobre numpy broadcasting.

4

En lugar de escribir una nueva función, numpy.ix_ debe hacer lo que quiera.

Aquí se muestra un ejemplo de la documentación:

>>> ixgrid = np.ix_([0,1], [2,4]) 
>>> ixgrid 
(array([[0], 
    [1]]), array([[2, 4]])) 
>>> ixgrid[0].shape, ixgrid[1].shape 
((2, 1), (1, 2))' 
+4

sería bueno si pudieras decir cómo? ... –

4

He aquí una versión multidimensional de meshgrid que escribí:

def ndmesh(*args): 
    args = map(np.asarray,args) 
    return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)]) 

Tenga en cuenta que las matrices devueltas son vistas de la matriz de datos originales, así que cambiar las matrices originales afectará a las matrices de coordenadas.

+0

hizo el truco para mí, muchas gracias! – somada141

31

Numpy (a partir de 1.8 creo) ahora es compatible con una mayor generación 2D de cuadrículas de posición con meshgrid. Una adición importante que realmente me ayudó es la capacidad de elegir el orden de indexación (ya sea xy o ij para la indexación cartesiana o matriz, respectivamente), lo que verifiqué con el siguiente ejemplo:

import numpy as np 

x_ = np.linspace(0., 1., 10) 
y_ = np.linspace(1., 2., 20) 
z_ = np.linspace(3., 4., 30) 

x, y, z = np.meshgrid(x_, y_, z_, indexing='ij') 

assert np.all(x[:,0,0] == x_) 
assert np.all(y[0,:,0] == y_) 
assert np.all(z[0,0,:] == z_) 
Cuestiones relacionadas