Fixing WCS bug

[ConnorStoneAstro]

@wmwv identified an issue in #115 which ultimately was tracked to the way AstroPhot converts world coordinates (RA,DEC) into the tangent plane where AstroPhot does it's calculations. It turns out that to a large degree the original coordinate system was implicitly defined and the conversion between world coordinates and tangent plane coordinates only really worked on the equator.

The new system makes explicit the position where the tangent plane and the celestial sphere (world coordinates) meet. The update also provides clear functions to map world to/from tangent plane and also tangent plane to/from world. A new WCS base class added for Window and Image_Header gives them the ability to call function that map world coordinates. The new class maintains a global reference (RA_0, DEC_0) so that multi-band multi-epoch images can be mapped to the same tangent plane easily.

This also comes with some nomenclature changes. The old world_to_pixel function(s) was pulling double duty and just causing problems so now it's world_to_plane and plane_to_pixel. Docstrings have been updated, though I may have missed some.

A docs page is added to go into more detail. This can be found at docs/coordinates.rst and when merged it will be a new page in the docs website. It may be beneficial to create a notebook tutorial on coordinate systems for users that care deeply about exact coordinates.

I have also updated the tutorial notebooks correspondingly at: https://github.com/Autostronomy/AstroPhot-tutorials/tree/wcsbug

This PR closes #115 and would constitute a version update breaking backwards compatibility. This will be version 0.12.0 unless I get the parameter DAG update done before we are happy with this WCS update.

[codecov[bot]]

Codecov Report

Merging #116 (d773fb1) into main (90ce051) will increase coverage by 1.68%. The diff coverage is 94.23%.

@@            Coverage Diff             @@
##             main     #116      +/-   ##
==========================================
+ Coverage   75.13%   76.81%   +1.68%     
==========================================
  Files          77       78       +1     
  Lines        6804     6953     +149     
==========================================
+ Hits         5112     5341     +229     
+ Misses       1692     1612      -80     

Files	Coverage Δ
astrophot/image/__init__.py	100.00% <100.00%> (ø)
astrophot/image/jacobian_image.py	76.05% <100.00%> (+0.71%)	:arrow_up:
astrophot/image/model_image.py	68.57% <100.00%> (+1.90%)	:arrow_up:
astrophot/image/psf_image.py	94.87% <100.00%> (+0.13%)	:arrow_up:
astrophot/image/target_image.py	69.80% <100.00%> (ø)
astrophot/models/_shared_methods.py	85.66% <100.00%> (ø)
astrophot/models/airy_psf.py	92.30% <100.00%> (ø)
astrophot/models/edgeon_model.py	94.04% <100.00%> (ø)
astrophot/models/galaxy_model_object.py	95.83% <100.00%> (ø)
astrophot/models/model_object.py	81.16% <100.00%> (ø)
... and 8 more

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[ConnorStoneAstro]

Looks like I need to add more unit tests, that will be my goal on Monday.

[ConnorStoneAstro]

I've been thinking about the world to tangent plane reference coordinate. Right now I make that a class variable for the WCS class, effectively making it a global variable. This is nice since it forces everything to be on the same tangent plane. But it is bad since anyone tinkering around can get surprised by the behavior which changes depending on which image is created first. So perhaps a better option would be to let each object have it's own reference coordinates, however if someone tries to construct an Image_List object it will throw an error if they are not all on the same tangent plane.

That seems more pythonic to me, so perhaps I'll refactor to make that happen. It's actually not a very big change so it shouldn't change much else.

[ConnorStoneAstro]

Ok, so now every image can have it's own reference RA, DEC. Looking closer at the fits standard it seems they keep track of the reference RA, DEC and also the point in the image that the reference corresponds to. I think I was doing this indirectly and only allowing the image center or bottom corner to be the reference point. I'll make that more general soon so there is no ambiguity.

[ConnorStoneAstro]

Converted to draft since it's clear it still needs refining and more tests.

[ConnorStoneAstro]

Ok, I have added enough generality for users to very clearly specify where images go in relation to each other. However, the Window objects now can't keep up with this level of generality. So next I'll need to make it so they can handle having a pivot point other than the window origin. After that, the new more comprehensive coordinate system should be ready to go!

[wmwv]

Good luck! I've been following along and checking things out along the way.

[ConnorStoneAstro]

Thanks for checking things out along the way :)

I've been trying to reconcile the Window objects with the new coordinate capabilities. My original vision for the Window objects was that they would purely represent coordinates on the tangent plane. However, now that I am using the pixelscale matrix, and reference points to define all the coordinate systems, it seems forced to keep it like this. I think instead I'll make windows more like Astropy WCS objects in that they define all the transformations from world coordinates to pixel coordinates. This means I can mostly just move capabilities over from Image_Header to Window and then everything should behave nicely (fingers crossed).

[wmwv]

Sounds good. I'm working on putting together a little test for something that's currently failing. (The failing is easy; making it the minimally reproducible test case is harder.)

[ConnorStoneAstro]

Sounds good. I'm working on putting together a little test for something that's currently failing. (The failing is easy; making it the minimally reproducible test case is harder.)

If the failing is something in how AstroPhot works, fell free to detail it here or on the issue and I can try to build the solution into this PR.

[wmwv]

Right, didn't mean to be mysterious, I just can't write a simple test case yet.

The issue is that the default window doesn't always actually match the pixel values of the image.

The function that generates the default window is not a separate function, but is rather wrapped into the image_header.init so I'm working on creating a minimal header that captures this.

[ConnorStoneAstro]

Ah, I see. Yes I think this is the issue I'm hoping to clear up with the next couple commits

[wmwv]

Here's an example of a test I added to tests/test_window.py to test a set of pixelscale and nx, ny that I get from an actual image that fails to ap.image.plot because the default window size is wrong. The default window size gets pixels of 4075, 3999 instead of 4072, 4000.

    def test_window_roundtrip(self):

        nx, ny = 4072, 4000
        pixelscale = torch.tensor(
            [[-0.1872, 0.0702], [0.0701, 0.1870]],
            dtype=torch.float64)
        origin = [288.3701, -531.1276]
        end = [-482.0643, 1033.4533]
        window1 = ap.image.Window(origin=origin, shape=end)
        indices = window1._get_indices(window1, pixelscale)

        self.assertTrue(indices[0].stop.item() == ny)
        self.assertTrue(indices[1].stop.item() == nx)

So, I don't exactly expect this test to pass in anything you fix. But rather I anticipate your fixes will result in a pixelscale, origin, and end that correctly recover nx, ny.

[wmwv]

Here's the WCS string for the relevant image. Sorry I can't format it better, but the FITS standard defines exact character chunks and I don't get things right when I try to manually line break to make the code readable.

wcs_header_string_nnnp = "WCSAXES =                    2 / Number of coordinate axes                      CRPIX1  =          2015.192126 / Pixel coordinate of reference point            CRPIX2  =          1968.704597 / Pixel coordinate of reference point            PC1_1   = -5.5410852767805E-05 / Coordinate transformation matrix element       PC1_2   = -7.1379810679143E-07 / Coordinate transformation matrix element       PC2_1   = -7.0957503843761E-07 / Coordinate transformation matrix element       PC2_2   =  5.5515960689837E-05 / Coordinate transformation matrix element       CDELT1  =                  1.0 / [deg] Coordinate increment at reference point  CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  CUNIT1  = 'deg'                / Units of coordinate increment and value        CUNIT2  = 'deg'                / Units of coordinate increment and value        CTYPE1  = 'RA---TAN'           / TAN (gnomonic) projection + SIP distortions    CTYPE2  = 'DEC--TAN'           / TAN (gnomonic) projection + SIP distortions    CRVAL1  =      60.371857275079 / [deg] Coordinate value at reference point      CRVAL2  =     -44.125292354859 / [deg] Coordinate value at reference point      LONPOLE =                180.0 / [deg] Native longitude of celestial pole       LATPOLE =     -44.125292354859 / [deg] Native latitude of celestial pole        MJDREF  =                  0.0 / [d] MJD of fiducial time                       RADESYS = 'ICRS'               / Equatorial coordinate system                   END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             "

[wmwv]

@ConnorStoneAstro Am I correct that you're basically on a path to re-implementing a WCS library in pytorch? Is it worth being explict about this and making this part a clearly separate module in the code? (not yet a separate package -- although one could imagine a future in which people find that useful).

[ConnorStoneAstro]

More or less, yes I am working in that direction. The scope of this PR kind of expanded from what I originally intended. I could see this becoming it's own separate module like you mention (maybe one day a package, but probably not). There is a major difference from Astropy WCS though in that between the world coordinates and pixel coordinates I have the tangent plane. My goal is for multi-image analysis that the tangent plane is where everything can be matched up between images. I'm not sure this is a feature in Astropy WCS since it goes directly between world and pixel coordinates

[ConnorStoneAstro]

Actually my latest plan with the Window objects absorbing the pixel coordinates will basically make them a stand alone thing that may be useful elsewhere.

[wmwv]

I think AstroPy WCS does provide the same intermediate representation through proj_plane_pixe_scales: https://docs.astropy.org/en/stable/api/astropy.wcs.utils.proj_plane_pixel_scales.html But I'm not 100% sure if it's how you're thinking about it.

After you finish this PR, you'll probably be well-primed to read this series of WCS papers by Greisen and Calabretta talking about defining convenient intermediate representations between the sky and a detector: https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1061G https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1077C

• Fig. 1 of each paper is particularly relevant in understanding the framework the authors are thinking about

And then when you've mastered that, you can go on to thinking about WCS as a way to encode spectral information :) https://ui.adsabs.harvard.edu/abs/2006A%26A...446..747G/abstract

[ConnorStoneAstro]

Ah, looking at those references, especially: https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1077C that is exactly where I am headed! Once again I find myself re-inventing the wheel. Oh well, I've learned a lot along the way and I ended up at the same place. This gives me some confidence though. And maybe I will try to separate the logic out a bit.

[ConnorStoneAstro]

Ok saving progress for now. This is a major refactor. I've collected the WCS stuff into a standalone WCS class. The interface with WCS is now almost entirely through window objects. I haven't worked out all the bugs yet. But I am going to a wedding this week so I might not get much more done. I've set aside all of next week to finish this and implement the parameter DAG.

It's getting close! Still some indexing bugs to iron out though. Looking back, there wasn't actually a "bug" per-se in AstroPhot before, I was just conflating world coordinates and tangent plane coordinates. But now I've come this far I really like the new system (at least once its working). It should be really easy to just give it a WCS object and go from there.

[ConnorStoneAstro]

Had some time to work on it. It now seems to be running correctly, though I've cut out a bunch of (I think) superfluous functionality. Going to add new tests that check what I actually want the WCS stuff to do.

At this point I think I've added all the needed WCS functionality to fairly seamlessly work with celestial coordinates and astropy WCS objects. Also, the WCS stuff has all been abstracted away into its own class so (hopefully) it wont be too painful later to add distortion parameters.

[wmwv]

Yeah, I started a little thinking about where SIP might fit in to the new WCS class. Sip is pixel<->focal plane<->projection plane<->sky(ra, dec) instead of the current pixel<->projection plane<->sky(ra, dec)

I don't know if this should be another option under the PPCS class or an added class.

[ConnorStoneAstro]

Ok, I think this is ready to go! Sorry the scope of this PR grew well beyond the initial problem and is now at over 2,000 affected lines... I'm pretty happy with the results though!

To your comment about adding SIP, I think the way to do it with the minimum change would be to make it an option for the PPCS. However, since everything is abstracted under the WCS class, we could restructure the hierarchy if making the SIP a distinct class seems worthwhile.

[ConnorStoneAstro]

Found an indexing bug, but didn't have time to fix it today. Changed back to draft

[wmwv]

I've started an issues for the WCS SIP #118 and started an issues/118 branch to work on it.

[ConnorStoneAstro]

I think this PR has enough in it already, but perhaps a future update for the built-in plotting methods could map to world coordinates then use a WCS object to plot images with RA/DEC projection and also include NE arrows :)

[ConnorStoneAstro]

Ok, unit tests run and also running the tutorial notebooks seems to work too. Now I think it's finally ready!

[ConnorStoneAstro]

Ok, I'll try to make all the adjustments and then merge the updates :) Thanks for going through everything!