Autostronomy / AstroPhot

A fast, flexible, automated, and differentiable astronomical image 2D forward modelling tool for precise parallel multi-wavelength photometry
https://astrophot.readthedocs.io
GNU General Public License v3.0
88 stars 9 forks source link

Creating a window with a certain WCS doesn't produce a positive shape. #115

Closed wmwv closed 1 year ago

wmwv commented 1 year ago

Describe the bug Using a window object based on a WCS doesn't work when constructing a Target_Image object

WCS """ WCS Keywords

Number of WCS axes: 2 CTYPE : 'RA---TAN-SIP' 'DEC--TAN-SIP'
CRVAL : 60.1432867246932 -44.1329261633904
CRPIX : 2053.100183 1971.914889
CD1_1 CD1_2 : -5.06818049270775e-05 2.26561792162743e-05
CD2_1 CD2_2 : 2.26423678359844e-05 5.05664537482149e-05
NAXIS : 4072 4000 """

target = ap.image.Target_Image(
        data=np.array(img, dtype=np.float64),
        variance=var,
        zeropoint=22.5,
        psf=psf,
        wcs=wcs,
    )

Creating window with

sn = {"ra": 60.2901401, "dec": -44.142051}
sn_coord = SkyCoord(sn["ra"], sn["dec"], unit=u.degree)
window_size = 5
window = ap.image.window_object.Window(
    center=(sn_coord.ra.arcsecond, sn_coord.dec.arcsecond),
    shape=(window_size, window_size),
    projection = target.pixelscale,
)

where target.pixelscale is

tensor([[-0.1825,  0.0816],
        [ 0.0815,  0.1820]], dtype=torch.float64)

To Reproduce Steps to reproduce the behavior:

  1. Go to '...'
  2. Click on '....'
  3. Scroll down to '....'
  4. See error

Expected behavior A clear and concise description of what you expected to happen.

Screenshots If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

Smartphone (please complete the following information):

Additional context Add any other context about the problem here.

wmwv commented 1 year ago
print(target_r_0.window.projection)
print(target_r_0.pixelscale / torch.linalg.det(target_r_0.pixelscale).abs().sqrt())

produces

tensor([[-0.9138,  0.4085],
        [ 0.4083,  0.9118]], dtype=torch.float64)
tensor([[-0.9138,  0.4085],
        [ 0.4083,  0.9118]], dtype=torch.float64)
ConnorStoneAstro commented 1 year ago

Thanks for catching this!

Ok, I think I reproduced the problem. It seems to be related to how out of boundary slices are handled. The window you requested at (5,5) arcsec and the center of the supernovae must have gone past the borders of the target image, possibly entirely outside the image. The proper behavior should be to just cut down the border to the image size in the appropriate dimensions, but it didn't do that for some reason. So now I just need to fix the behavior under that condition. Fingers crossed, it'll be fixed and pushed tonight.

Here I have reproduced the negative shape (and later it crashes trying to plot) Screenshot from 2023-09-15 17-41-33

Also added a repr, so that will be available soon too

ConnorStoneAstro commented 1 year ago

Observe that requesting a window fully inside the image doesn't have the error.

Screenshot from 2023-09-15 17-45-54

Screenshot from 2023-09-15 17-46-22

This is still very much a bug I need to fix! At the very least it should throw a helpful error message if you try to request a window that is outside the image.

wmwv commented 1 year ago

I’ll be happy to review the PR, if you like.

ConnorStoneAstro commented 1 year ago

That would be great!

Also, I'm sure you were giving coordinates correctly. I suspect that the reason your supernova coordinates were landing outside the image is because the sky is a sphere and your image was near enough to a pole where RA is modified by cos(dec). So what I really need is a way to communicate with image and window objects that explicitly uses RA and DEC and moves those into the assumed tangent plane. This means adding a few fairly simple functions. Fingers crossed sometime tonight!

ConnorStoneAstro commented 1 year ago

Ok, finishing tonight won't be happening, but I have all day tomorrow to work on this. I've been thinking a lot on how to do this "properly". I'm going to stick with a tangent plane, though perhaps this is something to revisit later.

I briefly considered actually doing everything in spherical coordinates, but this increases the number of calculations for every model since R = sqrt((x2-x1)^2 + (y2-y1)^2) becomes R = arccos(sin(x1)sin(x2) + cos(x1)cos(x2)cos(y2 - y1)) and apparently this formula is numerically unstable for small R which is the regime I'm mostly in. So that won't work.

Ok, so I'd like to do a tangent plane, but I have to deal with some problems. First, off the equator, RA is modified by cos(dec) so that 1 arcsec in RA is not 1 arcsec on the sky anymore. In the worst case you may have a celestial pole in your image and then there's a coordinate singularity... Second, in these cases off the equator it can be hard to get multiple images on the same tangent plane since their reference pixels are not necessarily at the same real coordinate, meaning that the P\vec{v} + \vec{o} (P is pixelscale matrix, v is position vector, o is origin) methodology will not work so well.

My solution is to have a parameter reference_radec which is explicitly the (RA_0, DEC_0) coordinates at which the tangent plane contacts the celestial sphere. That way a user provides coordinates in (RA, DEC) and I can do a Gnomonic projection onto the tangent plane. The (0,0) point of the tangent plane will be where it contacts the sphere, north (at the contact point) will be positive y-axis, west (at the contact point) will be positive x-axis. This way, multiple images can be projected onto the same tangent plane in a (mostly) consistent way. To make this explicit I think I'll go through everything and change pixel_to_world into pixel_to_plane and so on. This way the meaning of the functions will be clearer. And then I can maybe include pixel_to_world functions that will actually map to (RA, DEC) that users can take advantage of.

This message is a log for myself and an update on my thinking in case there is anything you think I've missed!

ConnorStoneAstro commented 1 year ago

Quick update. I have implemented all the projection functions for Gnomonic, Orthographic, and Steriographic projections depending on what the user prefers.

Here are a bunch of longitude lines projected in the various systems: large_field_longitude

My reference has been wolfram which describes each system in a consistent framework (gnomonic, orthographic, steriographic).

And here are some centered longitude lines of length 2 arcseconds projected by each system. They all are locally flat enough that no difference can be seen between them. small_field_longitude

Now I just need to build out the interface and add docstrings to make it clear what is going on internally.

wmwv commented 1 year ago

I'm impressed to see such progress in < 24 hours!

I agree with the tangent plane.

wmwv commented 1 year ago

🎉

ConnorStoneAstro commented 1 year ago

🎉 thanks @wmwv !