Better handling of 2D and 3D scenarios

[bwinkel]

The longer I use cygrid, the more inconvenient it feels that cygrid always needs a third axis, even for continuum data. This often means, one needs to create a copy of the (2D) FITS header, and add a third dimension, before it can be plugged to the WcsGrid constructor. Likewise, the input data array would need an explicit reshape to add the spectral axis, even though it is redundant.

The reasons for this design choice are that the underlying Cython/C++ code cannot easily handle numpy broadcasting rules. Therefore, my guess would be that under the hood it will always be simpler to assume the 3D case. However, one could perhaps change the higher-level handlers to add the 3rd dimension on-the-fly to relieve the burden from the end user.