title: How to extract a 2D plane from a 3D numpy meshgrid
body: <p><strong>[TLDR]:</strong></p>
  
  <p>Essentially my question boils down to how one can extract the 2d data of a  plane from a 3D numpy meshgrid</p>
  
  <p><strong>[Detailed Description]:</strong></p>
  
  <p>I am calculating the electric field of two (or more) point charges. I did this in 2D and can plot the results via matplotlib using quiver or streamplot</p>
  
  <pre><code>import numpy as np
  from matplotlib import pyplot as plt
  
  eps_0 = 8e-12
  fac = (1./(4*np.pi*eps_0))
  
  charges  = [1.0,-1.0]
  qx       = [-2.0,2.0]
  qy       = [0.0,0.0]
  
  # GRID
  gridsize = 4.0
  N = 11
  X,Y = np.meshgrid( np.linspace(-gridsize,gridsize,N),
                     np.linspace(-gridsize,gridsize,N))
  # CALC E-FIELD   
  sumEx = np.zeros_like(X)
  sumEy = np.zeros_like(Y)
  
  for q, qxi, qyi in zip(charges,qx,qy):
      dist_vec_x = X - qxi
      dist_vec_y = Y - qyi 
      dist = np.sqrt(dist_vec_x**2 + dist_vec_y**2)
  
      Ex = fac * q * (dist_vec_x/dist**3)
      Ey = fac * q * (dist_vec_y/dist**3)
  
      sumEx += Ex
      sumEy += Ey
  
  # PLOT
  fig = plt.figure()
  ax = fig.add_subplot(111)
  ax.streamplot(X,Y,sumEx,sumEy)
  plt.show()
  </code></pre>
  
  <p>This produces the correct results
  <a href="https://i.stack.imgur.com/gRo1q.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/gRo1q.png" alt="2D dipole"></a></p>
  
  <p>I can easily extend this to 3D</p>
  
  <pre><code>import numpy as np
  from mpl_toolkits.mplot3d import Axes3D
  from matplotlib import pyplot as plt
  
  eps_0 = 8e-12
  fac = (1./(4*np.pi*eps_0))
  
  charges = [1.0,-1.0]
  qx      = [-2.0,2.0]
  qy      = [0.0,0.0]
  qz      = [0.0,0.0]
  
  # GRID
  gridsize = 4.0
  N = 11
  X,Y,Z = np.meshgrid( np.linspace(-gridsize,gridsize,N),
                       np.linspace(-gridsize,gridsize,N),
                       np.linspace(-gridsize,gridsize,N))
  
  # CALC E-FIELD   
  sumEx = np.zeros_like(X)
  sumEy = np.zeros_like(Y)
  sumEz = np.zeros_like(Z)
  for q, qxi, qyi, qzi in zip(charges,qx,qy,qz):
      dist_vec_x = X - qxi
      dist_vec_y = Y - qyi
      dist_vec_z = Z - qzi
  
      dist = np.sqrt(dist_vec_x**2 + dist_vec_y**2 + dist_vec_z**2)
  
      Ex = fac * q * (dist_vec_x/dist**3)
      Ey = fac * q * (dist_vec_y/dist**3)
      Ez = fac * q * (dist_vec_z/dist**3)
  
      sumEx += Ex
      sumEy += Ey
      sumEz += Ez  
  
  # PLOT
  fig = plt.figure()
  ax = fig.gca(projection='3d')
  ax.quiver(X,Y,Z,sumEx,sumEy,sumEz, pivot='middle', normalize=True)
  plt.show()
  </code></pre>
  
  <p>This also yields the correct result when plotted in 3D (as far as I can tell)</p>
  
  <p><a href="https://i.stack.imgur.com/8goWm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8goWm.png" alt="3D dipole"></a></p>
  
  <p>But for some reason I can not figure out how to extract the data from one x-y plane from the generated 3D numpy mesh. I thought I could just do something like</p>
  
  <pre><code>zplane = round(N/2)
  ax.quiver(X,Y,sumEx[:,:,zplane],sumEy[:,:,zplane])
  </code></pre>
  
  <p>but this does not do the trick. Does anyone know the proper way here?</p>