Order of indices in 3D mesh tally

I have come across something a bit strange. It appears that there is a shift in the dimensions when trying to extract data from a 3D-mesh tally.

I define a mesh tally like so:

mesh_dim=(8,8,15)
meshxyz=openmc.RegularMesh()
meshxyz.lower_left=[-6,-6,-6]
meshxyz.upper_right=[6,6,6]
meshxyz.dimension=mesh_dim
mf=openmc.MeshFilter(meshxyz)
map_tally=openmc.Tally(name='flux')
map_tally.scores=['flux']

What I’d like to do (among other things) is to integrate the resulting tally along e.g. z.
I tried to do something akin to:

sp=openmc.StatePoint('statepoint.10.h5')
sl=sp.get_tally(name='flux').get_slice(['flux'])
sl.mean.shape=(8,8,15)
A=sl.mean.sum(axis=2)

But this matrix looked like something was off with the indices. So I tried to change the last two lines to:

sl.mean.shape=(8,15,8)
A=sl.mean.sum(axis=1)

… which gave me the expected results. The fact that the axis order appears to be different here seems a bit counterintuitive to me. Does anyone have any suggestions as to what is going on here? I could go to a dataframe formulation but I’d rather not.

As an add-on - I also added energy resolution to my tally so it became 4D. by adding these lines.

ef=openmc.EnergyFilter(np.logspace(3,7,20))
map_tally.filters=[mf,ef]

This on the other hand did what I had expected, as the energy axis became the 4th dimension.
I.e. to get back to the same plot as before I’d do:

sl.mean.shape=(8,15,8,19)
A=sl.mean.sum(axis=(1,3))