The current approach is correct. It follows ITK conventions ( and ITK Python ), and numpy conventions. ITK index via [i,j,k] and numpy index via [k,j,i]. Look at other Python imaging libraries such as PIL or scikit-image.
Yes you can utilize numpy.reshape to convert to fortran order while reversing the shape array, but that may result in a copy. And when you apply algorithms to the fortran array it may have to copy the data so that the buffer is in the proper order. This is not recommended.
Regarding keras, you have conflated index order and image data layout options. keras is following numpy conventions for [k,j,i] indexing. For RGB ( multi-components ), it has the option for the RGB to be the fast axis or the slow axis. For the former you can consider the layout as each pixel is a RGB value, while with the latter the data is laid out as the R-image, then G-image, then B-image.