can explore with a python ast interface

the AST user manual
8.8 Representing Offsets from a Specified Sky Position
provides some clues.

As an alternative (instead of changing the entire grid coordinate system), we can use the trick from the region projection by generating a big "cross" with consideration of projection effect as an overlay on the image.

Solution found (via AST directly)

To do so, we can make a fake header with CTYPE1/2 as "LINEAR", then the axes will become a linear one. The custom reference pixel can be set with CRPIX1/2. The delta angular distance is referred from CDELT1/2 as a constant. Therefore, if the raw image has projection distortion, the computed offset axes will have significant errors (because the axes are nonlinear actually). If we adopt this trick, we will need to pop a warning to users when the image has visible projection distortion.