What's the difference between the wcs package in astropy and LSST pipeline

Support
#1

Recently I meet some questions while dealing with ZTF data using your pipeline.

First is about wcs,
while calibrating, the terminal will output the following:

calibrate INFO: Loading reference objects from region bounded by [7.08789730, 7.67838167], [61.66964322, 61.94886698] RA Dec
calibrate INFO: Loaded 6320 reference objects

And it’s related to the loadReferenceObjects.py, what it do is to read the region of the image, and load a refCat region a little larger than that.
But while using the wcs package in astropy to inspect the distribution of the ra,dec of the image. It gets a different result, for example:

from astropy.io import fits
from astropy.wcs import WCS
hdu = fits.open("/home/yu/ztf/data/raw/after_split/ztf_20210804364132_c11_41_zr.fits")
wcs = WCS(hdu[1].header)
wcs.pixel_to_world_values(10,10)

and it output:

(array(4.52203392), array(63.95023313))

which is not in the region [7.08789730, 7.67838167], [61.66964322, 61.94886698], what caused this? Is there difference between the wcs package in astropy and LSST pipeline?

Thank you!

#3

One immediate thing that comes to mind is that sometimes, if a given image is a piece (a patch) of a larger area (a tract), the pixels do not necessarily start at 1,1. So if pixel 10,10 is not in the image, the WCS could be working fine but will return RA,Dec that also are not in the image. So that’s at least one first thing to check.

#4

Thank you for your help!

I guess it’s not because the image is a piece of a larger area and there are some detailes about the question:

Inspect what in collection:

for ref in registry.queryDatasets('raw',collections='WFST/raw/all',instrument='WFST'):
    print(ref.dataId.full)
{band: 'g', instrument: 'WFST', detector: 4, physical_filter: 'WFST-G', exposure: 1}
{band: 'g', instrument: 'WFST', detector: 4, physical_filter: 'WFST-G', exposure: 2}
{band: 'g', instrument: 'WFST', detector: 4, physical_filter: 'WFST-G', exposure: 3}
{band: 'g', instrument: 'WFST', detector: 4, physical_filter: 'WFST-G', exposure: 11}

the header of the raw is:

SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                 9216                                                  
NAXIS2  =                 9232                                                  
EXTEND  =                    T                                                  
CRPIX1  =               4608.5                                                  
CRPIX2  =               4616.5                                                  
CRVAL1  =                  180                                                  
CRVAL2  =                   36                                                  
CD1_1   = 9.16666666666666E-05                                                  
CD1_2   =                    0                                                  
CD2_1   =                    0                                                  
CD2_2   = 9.16666666666666E-05                                                  
CTYPE1  = 'RA---TAN'                                                            
CTYPE2  = 'DEC--TAN'                                                            
EXPTIME =                   30                                                  
DARKTIME=                   30                                                  
AIRMASS =                    1                                                  
EXP-ID  = 'WFST00000060'                                                        
OBJECT  = '9999    '                                                            
CCDNUM  =                    4                                                  
HUMIDITY=                   20                                                  
OUTTEMP =                   20                                                  
PRESSURE=                621.6                                                  
INSTRUME= 'WFST    '                                                            
FILTER  = 'WFST-G  '                                                            
DATE-OBS= '2022-09-27'                                                          
UT-STR  = '10:55:27.728'                                                        
UT-END  = '10:55:57.728'                                                        
DATA-TYP= 'object  '                                                            
RA2000  =                  180                                                  
DEC2000 =                   36                                                  
ALTITUDE=                    0                                                  
AZIMUTH =                    0

Use Butler to get one of the raw data and get a wcs:

raw = butler.get('raw',visit=3,instrument='WFST',detector=4,exposure=3,collections='WFST/raw/all')
expMd = ref_match._getExposureMetadata(raw)
wcs_lsst = expMd.wcs

and the wcs is(Here still meet the problem of the pixel scale, pixel scale is 0.175665 arcsec/pixel, not the wanted value 0.33(3600*9.16e-5), why?)

FITS standard SkyWcs:
Sky Origin: (180.0000000000, +36.0000000000)
Pixel Origin: (4607.5, 4615.5)
Pixel Scale: 0.175665 arcsec/pixel

get bbox:

bbox = raw.getBBox()
bbox2D = geom.Box2D(bbox)

bbox2D is:

Box2D(minimum=Point2D(-0.5, -0.5), dimensions=Extent2D(9216, 9232))

Then get corners:

corners = bbox2D.getCorners()

[Point2D(-0.5, -0.5),
 Point2D(9215.5, -0.5),
 Point2D(9215.5, 9231.5),
 Point2D(-0.5, 9231.5)]

See the ra, dec of the corners:

[wcs_lsst.pixelToSky(corner) for corner in corners]

[SpherePoint(180.27707609528898*geom.degrees, 35.77449199938334*geom.degrees),
 SpherePoint(179.72292390471102*geom.degrees, 35.77449199938334*geom.degrees),
 SpherePoint(179.72133698717087*geom.degrees, 36.22486718513159*geom.degrees),
 SpherePoint(180.27866301282913*geom.degrees, 36.22486718513159*geom.degrees)]

Use the wcs of astropy to do the same thing:

hdu = fits.open("/home/yu/python_by_yuyu/mock_wfst_task/raw/WFST00000030.fits")
wcs_astropy = WCS(hdu[0].header)
wcs_astropy.pixel_to_world_values(-0.5,-0.5)

the wcs astropy get:

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 180.0  36.0  
CRPIX : 4608.5  4616.5  
CD1_1 CD1_2  : 9.16666666666666e-05  0.0  
CD2_1 CD2_2  : 0.0  9.16666666666666e-05  
NAXIS : 9216  9232

the RA, DEC astropy get:

(array(179.4806856), array(35.5757607))

which is not the same of Point2D(-0.5, -0.5) :

SpherePoint(180.27707609528898*geom.degrees, 35.77449199938334*geom.degrees

They even change in different directions
So, what casued this?

Thank you!

#5

For raw data how are you constructing the afw WCS object when you wrote the raw formatter for the butler? Are you building a WCS from the camera geometry and boresight like we do for LSSTCam? Your raw data doesn’t seem to define any spatial distortions at all.

Astropy is going to use the FITS header directly and won’t understand anything about distortions.