WCS world_to_pixel dependent on SkyCoord distance

[mkolopanis]

Description

I have been using WCS to add celestial objects to radio astronomy observations and noticed discrepancy between the orthographic projected coordinates of solar system bodies depending on if the distance to the observer is included in their SkyCoord or not.

The AltAz frame azimuth and altitude information for a given body (e.g. the moon) is the same whether the distance is included or not, but when projected into the WCS frame the body will be in significantly different parts of the image.

Moreover, I believe there may be a bug in the coordinate transformation when distance is included because for an observation at night (0400 UTC) for the Owen's Valley Radio Observatory Long Wavelength Array (OVRO-LWA), the projected image coordinates for the Sun exist when the body distance is included which is unphysical as the Sun has set by this point in the night.

Below is the FITS header in string form, I have also uploaded the FITS image to this google drive to be more easily accessible.

Attached below is the MWE pulled from a notebook where I produced this calculation.

FITS Header information

SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    4 / number of data axes                            
NAXIS1  =                 4096 / length of data axis 1                          
NAXIS2  =                 4096 / length of data axis 2                          
NAXIS3  =                    1 / length of data axis 3                          
NAXIS4  =                    1 / length of data axis 4                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
BSCALE  =                   1.                                                  
BZERO   =                   0.                                                  
BUNIT   = 'JY/BEAM '           / Units are in Jansky per beam                   
BMAJ    =                   0.                                                  
BMIN    =                   0.                                                  
BPA     =                   0.                                                  
EQUINOX =                2000. / J2000                                          
LONPOLE =                 180.                                                  
BTYPE   = 'Intensity'                                                           
TELESCOP= 'OVRO-LWA'                                                            
OBSERVER= 'lwams.py'                                                            
OBJECT  = 'zenith  '                                                            
ORIGIN  = 'WSClean '           / W-stacking imager written by Andre Offringa    
CTYPE1  = 'RA---SIN'           / Right ascension angle cosine                   
CRPIX1  =                2049.                                                  
CRVAL1  =     141.834089116459                                                  
CDELT1  =             -0.03125                                                  
CUNIT1  = 'deg     '                                                            
CTYPE2  = 'DEC--SIN'           / Declination angle cosine                       
CRPIX2  =                2049.                                                  
CRVAL2  =     37.1579901237008                                                  
CDELT2  =              0.03125                                                  
CUNIT2  = 'deg     '                                                            
CTYPE3  = 'FREQ    '           / Central frequency                              
CRPIX3  =                   1.                                                  
CRVAL3  =      63917724.609375                                                  
CDELT3  =            45937500.                                                  
CUNIT3  = 'Hz      '                                                            
CTYPE4  = 'STOKES  '                                                            
CRPIX4  =                   1.                                                  
CRVAL4  =                   1.                                                  
CDELT4  =                   1.                                                  
CUNIT4  = '        '                                                            
SPECSYS = 'TOPOCENT'                                                            
DATE-OBS= '2024-04-11T04:02:25.0'  

Expected behavior

Position of the Sun on projected into a ortho-slant image should be (nan, nan) at the observation time 2024-04-11T04:02:25.000 UTC at the location GeodeticLocation(lon=<Longitude -118.28166669 deg>, lat=<Latitude 37.23977727 deg>, height=<Quantity 1183.48 m>). Instead using astropy.coordinates.get_body and projecting to wcs coordinates results in pixel values of the sun being up.

How to Reproduce

#!/usr/bin/env python
# coding: utf-8

from pathlib import Path

import numpy as np
from astropy.coordinates import AltAz, Angle, EarthLocation, SkyCoord, get_body
from astropy.io import fits
from astropy.time import Time
from astropy.wcs import WCS

location = EarthLocation.from_geodetic(
    lat=Angle("37.239777271d"),
    lon=Angle("-118.281666695d"),
    height=1183.48,
)
input_file = Path("20240411_040229_highband-I-dirty.fits")


with fits.open(input_file) as hdu_list:
    header = hdu_list[0].header
    wcs = WCS(header)
    obstime = Time(header["DATE-OBS"], format="isot", scale="utc")
wcs_sky = wcs.slice(np.s_[0, 0])
ovro_altaz = AltAz(obstime=obstime, location=location)


# Print some information about our location, the image, and WCS
print(ovro_altaz)
print(header)
print(wcs_sky)


# Get the sun, but also a sun which has no known distance to our observer
sun = get_body("sun", time=obstime, location=location)
new_sun = sun.transform_to(ovro_altaz)
new_sun = SkyCoord(alt=new_sun.alt, az=new_sun.az, frame=ovro_altaz)


#Sun location is equivalent in the GCRS frame but does not contain the distance information
print(sun)
print(new_sun.transform_to(sun.frame))



# Altitude and Azimuth are also equivalent
print(sun.transform_to(ovro_altaz))
print(new_sun)


#but the projected pixel locations differ
# the "complete" sun which includes distance puts the Sun onto the image, but it has already set
# the "new_sun" is off the image as expected
print(wcs_sky.world_to_pixel(sun))
print(wcs_sky.world_to_pixel(new_sun)

Versions

import astropy
try:
    astropy.system_info()
except AttributeError:
    import platform; print(platform.platform())
    import sys; print("Python", sys.version)
    import astropy; print("astropy", astropy.__version__)
    import numpy; print("Numpy", numpy.__version__)
    import erfa; print("pyerfa", erfa.__version__)
    try:
        import scipy
        print("Scipy", scipy.__version__)
    except ImportError:
        print("Scipy not installed")
    try:
        import matplotlib
        print("Matplotlib", matplotlib.__version__)
    except ImportError:
        print("Matplotlib not installed")
## -- End pasted text --
Linux-6.6.26-1-MANJARO-x86_64-with-glibc2.39
Python 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
astropy 6.0.1
Numpy 1.25.2
pyerfa 2.0.1.4
Scipy 1.11.1
Matplotlib 3.8.2

[mkolopanis]

Thanks for picking this up! I still have my notebook server running and can verify that new_sun.obstime attribute is equivalent to the obstime object. But I am sure you will be able to see everything when you have time to look into this.

[astrofrog]

I think this is unrelated to WCS - one also gets different results if converted to FK5 which is one of the steps of putting it in the WCS:

I need a coordinates expert to chime in, but I think what is happening is that FK5 (and ICRS etc) are barycentric coordinate frames, so if you convert a 3D SkyCoord (including finite distance) into barycentric coordinates you won't get the RA/Dec that you would see the object at from Earth, but from the barycenter. This is equivalent to RA/Dec you would see from Earth if the object was at an infinite distance. This is why I think the coordinates without distance work fine (because the barycentric/topocentric difference doesn't matter then since distance is assumed to be infinite) and coordinates with distance don't.

[StuartLittlefair]

Your explanation is correct @astrofrog . One way of thinking about this is that for an observer at the barycenter the Sun would indeed be in the image at pixel (440, 1323).

I’m not sure this is independent of WCS though. This behaviour is surprising enough that there should be a warning if coordinates in a non-barycentric frame of reference is passed to world_to_pixel.

This is reminiscent of the feature that caused so much confusion is SkyCoord.separation when coordinate systems had different origins. Perhaps WCS should take advantage the machinery introduced in #16246 ?

[mkolopanis]

Thanks for looking into this. I'm a little confused. Am I interpreting what is being said correctly? When transforming a SkyCoord into a WCS, it always transforms through a barycentric coordinate frame and as a result, all locations in WCS are apparent ra/dec relative the sun-earth barycenter?

[astrofrog]

@mkolopanis - no it won't always transform to barycentric, but your WCS is in FK5 which is a barycentric system, and RA/Dec will be relative to the Sun/Earth barycenter (which is how RA and Dec are defined in FK5) - you see the same effect if you just convert your positions to RA/Dec without even a WCS (see #16430 (comment))

[mkolopanis]

Thanks for the reply, appreciate it. Yeah that makes a lot of sense. I did some more digging into the fits header after reading through all your messages earlier. Made me think if wsclean wrote out in another frame it out be different. Thanks for confirming.

…

On Fri, May 10, 2024, 16:53 Thomas Robitaille ***@***.***> wrote: @mkolopanis <https://github.com/mkolopanis> - no it won't always transform to barycentric, but your WCS is in FK5 which *is* a barycentric system, and RA/Dec will be relative to the Sun/Earth barycenter (which is how RA and Dec are defined in FK5) - you see the same effect if you just convert your positions to RA/Dec without even a WCS (see #16430 (comment) <#16430 (comment)>) — Reply to this email directly, view it on GitHub <#16430 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACB7BLO3YXTULRJHD2FOBWTZBVMXPAVCNFSM6AAAAABHPI3N32VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMBVGM4TANJRGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[github-actions]

Hi humans 👋 - this issue was labeled as Close? approximately 11 hours ago. If you think this issue should not be closed, a maintainer should remove the Close? label - otherwise, I will close this issue in 7 days.

If you believe I commented on this issue incorrectly, please report this here

[github-actions]

I'm going to close this issue as per my previous message, but if you feel that this issue should stay open, then feel free to re-open and remove the Close? label.

If this is the first time I am commenting on this issue, or if you believe I closed this issue incorrectly, please report this here