Add the ability to calculate the Jacaboian matrix to `astropy.wcs`.

### Description

A useful feature for handling coordinate transformations is the Jacobian matrix. Primarily, it allows for a local determination of how the transformation modifies area elements (by the eigenvalues/determinant) and basis orientation (by the eigenvectors). The prototypical application is to determine in the xy coordinates which directions are North and East. at a particular point in the image. In the simple case, this would just be reporting the elements of `CDi_j` correctly packed in an array, but when distortion is relevant the local Jacobian should be modified by this. This should correctly handle all of the special cases, though, like when the header specifies two scales and an angle.

Ideally, the signatures of the functions should follow the coordinate transform functions. Something like `Jac_all_world2pix`, with the same arguments and meanings as the `all_world2pix` function, but which returns an `N` by `NAXIS` by `NAXIS` array, instead of just `N` by `NAXIS`.

### Additional context

My primary use case is overlaying text on a ds9 image and getting the text to display upright. When the fits file does unusual things with the data axis order (i.e. the Jacobian has a reflection) and/or a 90 degree rotation, it can make the text display upside-down (even when the text angle is specified relative to the fk5 coordinates as 0.0).

It would also make it more straightforward to generate compass roses to overlay on non-ds9 images.

Start by reading the astropy.wcs coordinate-transform entry point all_world2pix and compare its arguments and meanings with the proposed Jac_all_world2pix API. Determine how the returned N by NAXIS by NAXIS array should represent CDi_j, distortions, scales, angles, axis order, and reflections, then validate the behavior for the stated text-overlay and compass-rose use cases.

Written by the indexing model from the issue text.