Incorrect reshape of vtk->numpy array?

Question

I am reading VTK uniform grid into python. When I visualize a slice through the data in Paraview, I get the following (correct) image:

Then I visualize the slice using via numpy & pylab using the following script:

import vtk
from vtk.util.numpy_support import vtk_to_numpy
import pylab

imr=vtk.vtkXMLImageDataReader()
imr.SetFileName('flow.vti')
imr.Update()
im=imr.GetOutput()

nx,ny,nz=im.GetDimensions()
orig=im.GetOrigin()
extent=im.GetExtent()
spacing=im.GetSpacing()

flowVtk=im.GetPointData().GetArray("|flow|")
flow=vtk_to_numpy(flowVtk).reshape(nx,ny,nz)
# bottom z-slice
flowZ0=flow[:,:,0]
# set extent so that axes units are physical
img=pylab.imshow(flowZ0,extent=[orig[0],orig[0]+extent[1]*spacing[0],orig[1],orig[1]+extent[3]*spacing[1]],cmap=pylab.gray())
img.set_clim(vmin=0,vmax=1000)
pylab.show()

which seems to be out-of-phase. I tried reordering dimensions in reshape(...), it did something, but it has never shown the data it is actually supposed to show.

Is there something obviously wrong?

EDIT: I also tried reshape((nx,ny,nz),order="F") (fortran ordering) and now I get a much better image (with jet colormap for better clarity) which is almost correct, but the data is suspiciously rotated by 90°, plus I would like some authoritative explanation which ordering to use and why (which one is used by VTK internally?).

EDIT2: to get the same view as in Paraview, I had to do pylab.imshow(np.rot90(flowZ0)); not sure why, so the question is still open: