WCS Information - Githubissues

kadrlica commented 3 years ago

It would be nice to carry around the WCS information from the original FITS files. This would allow you to reconstruct the WCS if/when you create individual images. This will be useful if/when we want to do any follow-up or additional investigation of any individual cutouts.

I've written a script that demonstrates how this information can be accessed.

#!/usr/bin/env python
__author__ = "Alex Drlica-Wagner"

import numpy as np
import pylab as plt
import fitsio
import pandas as pd

DTYPE = [ (f"CTYPE{i}",'S8') for i in [1,2]] +\
        [ (f"DCRPIX{i}",float) for i in [1,2]] + \
        [ (f"CRPIX{i}", float) for i in [1,2]]  + \
        [ (f"CRVAL{i}", float) for i in [1,2]]  + \
        [ (f"CD{j}_{i}",float) for i in [1,2] for j in [1,2]]   + \
        [ (f"PV{j}_{i}",float) for i in range(11) for j in [1,2]]

DTYPES = [ (n,(d,4)) for n,d in DTYPE]

def make_wcs_table(size,dtype=DTYPE):
    row =  np.zeros(size,dtype=dtype)
    row.fill(np.nan)
    return row

def get_wcs(filename,ext='SCI'):
    row = make_wcs_table(1)

    if not filename or str(filename) == 'nan': return row 

    hdr = fitsio.read_header(filename,ext=ext)
    for name,dt in DTYPE:
        try:
            row[name] = hdr[name] 
        except KeyError:
            pass

    return row

def get_all_wcs(filenames):
    rows = []
    for filename in filenames:
        rows.append(get_wcs(filename))

    return np.concatenate(rows)

dirname = 'schechter_cutouts_20210114'

f = fitsio.FITS(dirname+'/cutouts.fits')
meta = pd.read_csv(dirname+'/metadata.csv')

wcs_g = get_all_wcs(meta['FILEPATH_IMAGE_G'])
wcs_r = get_all_wcs(meta['FILEPATH_IMAGE_R'])
wcs_i = get_all_wcs(meta['FILEPATH_IMAGE_Z'])
wcs_z = get_all_wcs(meta['FILEPATH_IMAGE_I'])

wcs = make_wcs_table(len(meta),dtype=DTYPES)
for i,_wcs in enumerate([wcs_g,wcs_r,wcs_i,wcs_z]):
    for name,dt in DTYPE:
        wcs[name][:,i] = _wcs[name]

outfile = 'wcs.fits'
out = fitsio.write(outfile,wcs,clobber=True,extname='WCS')

print(fitsio.FITS(outfile))

rmorgan10 commented 3 years ago

Will try to implement this next week 👍

rmorgan10 commented 3 years ago

Something I'm noticing is I don't see DCRPIX1, DCRPIX2, or any of the PV names in the headers of the DES tiles, so pretty much everything raises a KeyError and just stays as a NaN value:

(des20a) [rmorgan@des81 Y6_Bulk_Coadd_Cutouts]$ python
Python 3.7.7 (default, Mar 26 2020, 15:48:22) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import fitsio
>>> fitsio.read_header('/data/des40.b/data/des/y6a2/coadd/image/DES2228+0001/DES2228+0001_r4575p01_i.fits.fz', ext='SCI')

XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / 8-bit bytes
NAXIS   =                    2 / 2-dimensional binary table
NAXIS1  =                   24 / width of table in bytes
NAXIS2  =                10000 / number of rows in table
PCOUNT  =             83819913 / size of special data area
GCOUNT  =                    1 / one data group (required keyword)
TFIELDS =                    3 / number of fields in each row
TTYPE1  = 'COMPRESSED_DATA'    / label for field   1
TFORM1  = '1PB(8557)'          / data format of field: variable length array
TTYPE2  = 'ZSCALE'             / label for field   2
TFORM2  = '1D'                 / data format of field: 8-byte DOUBLE
TTYPE3  = 'ZZERO'              / label for field   3
TFORM3  = '1D'                 / data format of field: 8-byte DOUBLE
ZIMAGE  =                    T / extension contains compressed image
ZTILE1  =                10000 / size of tiles to be compressed
ZTILE2  =                    1 / size of tiles to be compressed
ZCMPTYPE= 'RICE_ONE'           / compression algorithm
ZNAME1  = 'BLOCKSIZE'          / compression block size
ZVAL1   =                   32 / pixels per block
ZNAME2  = 'BYTEPIX'            / bytes per pixel (1, 2, 4, or 8)
ZVAL2   =                    4 / bytes per pixel (1, 2, 4, or 8)
ZSIMPLE =                    T / file does conform to FITS standard
ZBITPIX =                  -32 / number of bits per data pixel
ZNAXIS  =                    2 / number of data axes
ZNAXIS1 =                10000 / length of data axis 1
ZNAXIS2 =                10000 / length of data axis 2
ZEXTEND =                    T / FITS dataset may contain extensions
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
EXTNAME = 'SCI'                / 
EQUINOX =               2000.0 / Mean equinox
MJD-OBS =       56548.15423801 / Modified Julian date at start
RADESYS = 'ICRS'               / Astrometric system
CTYPE1  = 'RA---TAN'           / WCS projection type for this axis
CUNIT1  = 'deg'                / Axis unit
CRVAL1  =           337.108336 / World coordinate on this axis
CRPIX1  =               5000.5 / Reference pixel on this axis
CD1_1   =  -7.305555555556e-05 / Linear projection matrix
CD1_2   =                  0.0 / Linear projection matrix
CTYPE2  = 'DEC--TAN'           / WCS projection type for this axis
CUNIT2  = 'deg'                / Axis unit
CRVAL2  =             0.016667 / World coordinate on this axis
CRPIX2  =               5000.5 / Reference pixel on this axis
CD2_1   =                  0.0 / Linear projection matrix
CD2_2   =   7.305555555556e-05 / Linear projection matrix
EXPTIME =               1080.0 / Maximum equivalent exposure time (s)
GAIN    =       53.94172363784 / Maximum equivalent gain (e-/ADU)
SATURATE=       26321.58823898 / Saturation Level (ADU)
SOFTNAME= 'SWarp'              / The software that processed those data
SOFTVERS= '2.40.1'             / Version of the software
SOFTDATE= '2019-02-05'         / Release date of the software
SOFTAUTH= '2010-2012 IAP/CNRS/UPMC' / Maintainer of the software
SOFTINST= 'IAP  http://www.iap.fr' / Institute
AUTHOR  = 'unknown'            / Who ran the software
ORIGIN  = 'ccc0029.campuscluster.illinois.edu' / Where it was done
DATE    = '2019-06-20T22:32:44' / When it was started (GMT)
COMBINET= 'WEIGHTED'           / COMBINE_TYPE config parameter for SWarp
BUNIT   = 'COUNTS/S'           / units of pixels
FILTER  = 'i DECam SDSS c0003 7835.0 1470.0' / Unique filter identifier
BAND    = 'i'                  / Short name for filter
TILENAME= 'DES2228+0001'       / DES Tilename
TILEID  =                90149 / Tile ID for DES Tilename
RESAMPT1= 'LANCZOS3'           / RESAMPLING_TYPE config parameter
CENTERT1= 'MANUAL'             / CENTER_TYPE config parameter
PSCALET1= 'MANUAL'             / PIXELSCALE_TYPE config parameter
RESAMPT2= 'LANCZOS3'           / RESAMPLING_TYPE config parameter
CENTERT2= 'MANUAL'             / CENTER_TYPE config parameter
PSCALET2= 'MANUAL'             / PIXELSCALE_TYPE config parameter
DESFNAME= 'DES2228+0001_r4575p01_i.fits' / DES production filename
PIPELINE= 'multiepoch'         / DESDM pipeline name
UNITNAME= 'DES2228+0001'       / DESDM processing unit name
ATTNUM  =                    1 / DESDM processing attempt number
EUPSPROD= 'MEPipeline'         / eups pipeline meta-package name
EUPSVER = 'Y6A1+1'             / eups pipeline meta-package version
REQNUM  =                 4575 / DESDM processing request number
DES_EXT = 'IMAGE'              / DES name of the extension
FZALGOR = 'RICE_1'             / Compression type
FZDTHRSD= 'CHECKSUM'           / Dithering seed value
FZQVALUE=                   16 / Compression quantization factor
FZQMETHD= 'SUBTRACTIVE_DITHER_2' / Compression quantization method
RA_CENT =     337.108372527779 / RA center
DEC_CENT=    0.016630472222219 / DEC center
RAC1    =     337.473571640963 / RA corner 1
DECC1   =   -0.348562220896466 / DEC corner 1
RAC2    =     336.743100359037 / RA corner 2
DECC2   =   -0.348562220896466 / DEC corner 2
RAC3    =     336.743099004521 / RA corner 3
DECC3   =    0.381895543631536 / DEC corner 3
RAC4    =     337.473572995479 / RA corner 4
DECC4   =    0.381895543631536 / DEC corner 4
RACMIN  =     336.743099004521 / Minimum extent of image in RA
RACMAX  =     337.473572995479 / Maximum extent of image in RA
DECCMIN =   -0.348562220896466 / Minimum extent of image in Declination
DECCMAX =    0.381895543631536 / Maximum extent of image in Declination
CROSSRA0= 'N'                  / Does Image Span RA 0h (Y/N)
MAGZERO =                 30.0 / Mag Zero-point in magnitudes/s
HISTORY   Image was compressed by CFITSIO using scaled integer quantization:
ZQUANTIZ= 'SUBTRACTIVE_DITHER_2' / Pixel Quantization Algorithm
HISTORY     q = 16.000000 / quantized level scaling parameter
HISTORY   'SUBTRACTIVE_DITHER_2' / Pixel Quantization Algorithm
HISTORY   Thu Jun 20 18:24:12 2019 column_interp over mask 0x0001
ZDITHER0=                 3739 / dithering offset when quantizing floats
CHECKSUM= 'n1aGn0VEn0ZEn0ZE'   / HDU checksum updated 2019-06-21T04:24:29
DATASUM = '4243202715'         / data unit checksum updated 2019-06-21T04:24:29
>>>

rmorgan10 commented 3 years ago

The code above produces this if I do:

for name, dt in DTYPE:
    print(name, wcs[name])

CTYPE1 [[b'RA---TAN' b'RA---TAN' b'RA---TAN' b'RA---TAN']]
CTYPE2 [[b'DEC--TAN' b'DEC--TAN' b'DEC--TAN' b'DEC--TAN']]
DCRPIX1 [[nan nan nan nan]]
DCRPIX2 [[nan nan nan nan]]
CRPIX1 [[5000.5 5000.5 5000.5 5000.5]]
CRPIX2 [[5000.5 5000.5 5000.5 5000.5]]
CRVAL1 [[337.108336 337.108336 337.108336 337.108336]]
CRVAL2 [[0.016667 0.016667 0.016667 0.016667]]
CD1_1 [[-7.30555556e-05 -7.30555556e-05 -7.30555556e-05 -7.30555556e-05]]
CD2_1 [[0. 0. 0. 0.]]
CD1_2 [[0. 0. 0. 0.]]
CD2_2 [[7.30555556e-05 7.30555556e-05 7.30555556e-05 7.30555556e-05]]
PV1_0 [[nan nan nan nan]]
PV2_0 [[nan nan nan nan]]
PV1_1 [[nan nan nan nan]]
PV2_1 [[nan nan nan nan]]
PV1_2 [[nan nan nan nan]]
PV2_2 [[nan nan nan nan]]
PV1_3 [[nan nan nan nan]]
PV2_3 [[nan nan nan nan]]
PV1_4 [[nan nan nan nan]]
PV2_4 [[nan nan nan nan]]
PV1_5 [[nan nan nan nan]]
PV2_5 [[nan nan nan nan]]
PV1_6 [[nan nan nan nan]]
PV2_6 [[nan nan nan nan]]
PV1_7 [[nan nan nan nan]]
PV2_7 [[nan nan nan nan]]
PV1_8 [[nan nan nan nan]]
PV2_8 [[nan nan nan nan]]
PV1_9 [[nan nan nan nan]]
PV2_9 [[nan nan nan nan]]
PV1_10 [[nan nan nan nan]]
PV2_10 [[nan nan nan nan]]

kadrlica commented 3 years ago

DCRPIX1, DCRPIX2 are the DELTA offset in the CRPIX that comes from creating the cutout (i.e., it is the offset that needs to be applied to shift the WCS to the origin of the cutout image). I'm applying this offset directly to CRPIX[1,2] in stack2image, and we may consider doing that here as well. In that case, we don't need to carry around DCRPIX, though it might be helpful for provenance.

kadrlica commented 3 years ago

The PV terms are probably not needed for the TAN projection used by the coadd tiles. Do you think we should still carry them for the DES cutouts?

rmorgan10 commented 3 years ago

(sorry for the delay, just got back from vacation)

If you don't think the PV terms are needed then I think we should ignore them. I think as long as the origin is updated then that should be sufficient for the WCS

kadrlica commented 3 years ago

Sorry, somehow gmail decided github issue emails were spam...

The number of PV terms (and the form of the WCS generally) depends on the transformation applied to the image. I think that we should assume that we are only storing information that was in the original header. This will mean that cutouts derived from CCDs vs COADDS will have a different form for the WCS table, but it will make it programmatically simple to just say: "I'll take everything in this table and put it back in the header"

DES-SL / Y6_Bulk_Coadd_Cutouts

WCS Information #25