ejeschke / ginga

The Ginga astronomical FITS file viewer
BSD 3-Clause "New" or "Revised" License
121 stars 77 forks source link

wcs_astropy_ape14: Fall back to old way when no convergence #967

Closed pllim closed 2 years ago

pllim commented 3 years ago

Fix #960 . I don't understand why the "old way" works but not the APE 14 interface. I wonder if @mcara can advise here. The image is jbt7a3k7q_flc.fits.

Filename: jbt7a3k7q_flc.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     281   ()      
  1  SCI           1 ImageHDU       249   (4096, 2048)   float32   
  2  ERR           1 ImageHDU        56   (4096, 2048)   float32   
  3  DQ            1 ImageHDU        48   (4096, 2048)   int16   
  4  SCI           2 ImageHDU       245   (4096, 2048)   float32   
  5  ERR           2 ImageHDU        56   (4096, 2048)   float32   
  6  DQ            2 ImageHDU        48   (4096, 2048)   int16   
  7  HDRLET        1 NonstandardExtHDU     22   (60480,)      
  8  HDRLET        2 NonstandardExtHDU     22   (60480,)      
  9  HDRLET        3 NonstandardExtHDU     26   (112320,)      
 10  HDRLET        4 NonstandardExtHDU     26   (112320,)      
 11  HDRLET        5 NonstandardExtHDU     26   (112320,)      
 12  WCSCORR       1 BinTableHDU     59   14R x 24C   [40A, I, A, 24A, 24A, 24A, 24A, D, D, D, D, D, D, D, D, 24A, 24A, D, D, D, D, J, 40A, 128A]   
 13  WCSDVARR      1 ImageHDU        15   (64, 32)   float32   
 14  WCSDVARR      2 ImageHDU        15   (64, 32)   float32   
 15  D2IMARR       1 ImageHDU        15   (64, 32)   float32   
 16  D2IMARR       2 ImageHDU        15   (64, 32)   float32   
 17  WCSDVARR      3 ImageHDU        15   (64, 32)   float32   
 18  WCSDVARR      4 ImageHDU        15   (64, 32)   float32   
 19  D2IMARR       3 ImageHDU        15   (64, 32)   float32   
 20  D2IMARR       4 ImageHDU        15   (64, 32)   float32   
 21  HDRLET        6 NonstandardExtHDU     26   (112320,)      
 22  HDRLET        7 NonstandardExtHDU     26   (112320,)      
 23  HDRLET        8 NonstandardExtHDU     26   (112320,)  
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  
CRVAL : 168.9840199589  1.4973412658182  
CRPIX : 2048.0  1024.0  
CD1_1 CD1_2  : -2.1416145107558e-06  -1.3918446407705e-05  
CD2_1 CD2_2  : -1.3680995986142e-05  1.2042976985148e-06  
NAXIS : 4096  2048

This PR is a pre-requisite for #959 .

mcara commented 3 years ago

Fix #960 . I don't understand why the "old way" works but not the APE 14 interface. I wonder if @mcara can advise here.

On line 137 you call pix = self.wcs.wcs_world2pix(np.array([args], np.float_), origin) - notice wcs_world2pix - which does NOT use iteration to revert coordinates since it ignores polynomial (SIP) distortions and inversion can be performed analytically. APE14 calls all_world2pix which will fail if you are away by 1 degree from the center of an ACS image.

This may be an unnecessary complication, but instead of this complicated pirouette with adding 1 degree and catching exceptions and mixing APE14 with non-APE14 calls, could you create a bounding box around image(s) (multiple images for multiple ACS chips) and then create a spherical polygon and test if a point is inside the polygon before trying to invert coordinates? Just a suggestion.

mcara commented 3 years ago

Ah, adding +1 degree is in a different module - see https://github.com/ejeschke/ginga/issues/960#issue-909792583

I think that code could create spherical polygon and test if a point is inside (but I am not familiar with the code or why it is done the way it is done and so I may be wrong). Here it is OK to catch exceptions but, if you do not need logging all these messages, you could simply call all_world2pix(..., quiet=True) but you should test this to see if it works for you.

pllim commented 3 years ago

adding +1 degree is in a different module

I did try to make the "hack" smaller but it didn't fix the problem. Somewhere else in the code, it tries to do something that also triggers NoConvergence. Between Compass and the lower level wcs.py in the code base, I got lost where this is really happening. Maybe @ejeschke knows.

create spherical polygon

I hope you are not asking me to add spherical-geometry as a dependency here.

mcara commented 3 years ago

I hope you are not asking me to add spherical-geometry as a dependency here.

He-he. Yes, that's what I meant.

I did try to make the "hack" smaller but it didn't fix the problem.

What is the problem you are trying to solve by adding +1 degree?

mcara commented 3 years ago

Instead of spherical_geometry, you could take corners of the image and compute maximum "distance" (in degrees) from CRVAL and consider that to be the "radius of convergence" for coordinate inversion. If a point is inside the circle - call APE14, if not - call wcs_world2pix. Probably you can use something from astropy.coordinates for this.