Currently much of the functionality of this package, especially for the mesh manipulation, rely on an assumed matrix index ordering for spatial dimension and time. To make matters worse, that ordering is not strictly enforced across the package, and sometimes is projected onto the user to know what is what.

The xarray package (https://xarray.pydata.org) provides a nice, high-level abstract to allow dimensions to be more elegantly referenced by descriptive labels, which would greatly improve readability of the code, and reduce opportunity for user errors.

This is a target of future refactoring, but no core functionality is foreseen to change with this infrastructure (except for some API changes for some functions).