spacetelescope / STScI-STIPS

STScI-STIPS
https://stips.readthedocs.io
14 stars 16 forks source link

No PC or CD matrix in header #90

Closed benw1 closed 4 years ago

benw1 commented 4 years ago

Using version 1.0.8, there are no CD or PC matrix header keywords in the resulting images. WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 2048.0 / Pixel coordinate of reference point
CRPIX2 = 2048.0 / Pixel coordinate of reference point
CDELT1 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point
CDELT2 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection
CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection
CRVAL1 = 25.65 / [deg] Coordinate value at reference point
CRVAL2 = -37.21 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = -37.21 / [deg] Native latitude of celestial pole
DATEREF = '1858-11-17' / ISO-8601 fiducial time
MJDREFI = 0.0 / [d] MJD of fiducial time, integer part

york-stsci commented 4 years ago

Can you please provide the code that produced this image, and attach the produced FITS file. I'm unable to reproduce that issue with any of my local testing.

benw1 commented 4 years ago

Sorry for the delay. Here is the code. STIPS should also include, at least placeholders for the the full distortion information. Once STIPS can generate multi-chip images, we will need these values as well as the full distortion matrix (via 4th or 5th order SIP corrections) populated. For now, distortion is assumed to be zero. seed = np.random.randint(9999)+1000 scene_general = {'ra': 149.5, 'dec': 69.8, 'pa': pa, 'seed': seed} obs = {'instrument': 'WFI', 'filters': ['F062'], 'detectors': 1, 'distortion': False, 'oversample': 3, 'pupil_mask': '', 'background': 'avg', 'observations_id': dp_id, 'exptime': 1000, 'offsets': [{'offset_id': 4356, 'offset_centre': False, 'offset_ra': 0.0, 'offset_dec': 0.0, 'offset_pa': 0.0}]} obm = ObservationModule(obs, scene_general=scene_general, psf_grid_size=5) obm.nextObservation()

york-stsci commented 4 years ago

What are you using as the pa? The only thing I can think of is that astropy does not appear to generate a PC matrix if the x and y axes are exactly RA and DEC.

york-stsci commented 4 years ago

Okay, below is the code that I used (I made adjustments because the code provided above won't run because the variables 'pa' and 'dp_id' are not provided. In addition, when running this code I discovered a couple of bugs in the way that PSF caching interacts (or doesn't) with the config file, although I don't think that those issues could have affected the problem here.

from stips.observation_module import ObservationModule import numpy as np seed = np.random.randint(9999)+1000 pa = np.random.random()*360. scene_general = {'ra': 149.5, 'dec': 69.8, 'pa': pa, 'seed': seed} dp_id = np.random.randint(9999) obs = { 'instrument': 'WFI', 'filters': ['F062'], 'detectors': 1, 'distortion': False, 'oversample': 3, 'pupil_mask': '', 'background': 'avg', 'observations_id': dp_id, 'exptime': 1000, 'offsets': [ { 'offset_id': 4356, 'offset_centre': False, 'offset_ra': 0.0, 'offset_dec': 0.0, 'offset_pa': 0.0 } ] } obm = ObservationModule(obs, scene_general=scene_general, psf_grid_size=5) obm.nextObservation()

As far as I can see, this code doesn't produce any output files, so I certainly can't verify that there's no WCS information in whatever file(s) are supposed to be produced. Am I missing something?

benw1 commented 4 years ago

Sorry, it requires a few more lines: source_count_catalogues = obm.addCatalogue("catalog.txt") psf_file = obm.addError() fits_file, mosaic_file, params = obm.finalize(mosaic=False)

This requires a catalog file, which I attach here. You can set the dp_id and pa to 0. catalog.txt

york-stsci commented 4 years ago

When I run this, I get the following output:

>>> from astropy.io import fits >>> f = fits.open("sim_5596_000.fits") >>> f.info() Filename: sim_5596_000.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 4 () 1 SCA01 1 ImageHDU 55 (4096, 4096) float32 >>> f[0].header SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T >>> f[1].header XTENSION= 'IMAGE ' / Image extension BITPIX = -32 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 4096 NAXIS2 = 4096 PCOUNT = 0 / number of parameters GCOUNT = 1 / number of groups WCSAXES = 2 / Number of coordinate axes CRPIX1 = 2048.0 / Pixel coordinate of reference point CRPIX2 = 2048.0 / Pixel coordinate of reference point PC1_1 = 0.47570813752257 / Coordinate transformation matrix element PC1_2 = -0.87960318774707 / Coordinate transformation matrix element PC2_1 = 0.87960318774707 / Coordinate transformation matrix element PC2_2 = 0.47570813752257 / Coordinate transformation matrix element CDELT1 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point CDELT2 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection CRVAL1 = 149.5 / [deg] Coordinate value at reference point CRVAL2 = 69.8 / [deg] Coordinate value at reference point LONPOLE = 180.0 / [deg] Native longitude of celestial pole LATPOLE = 69.8 / [deg] Native latitude of celestial pole DATEREF = '1858-11-17' / ISO-8601 fiducial time MJDREFI = 0.0 / [d] MJD of fiducial time, integer part MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part RADESYS = 'ICRS' / Equatorial coordinate system EXTNAME = 'SCA01 ' / extension name DETECTOR= 'SCA01 ' FILTER = 'F158 ' DATE-OBS= '2020-09-30' TIME-OBS= 'Sep:09:1601488245' EQUINOX = 2000.0 PA_APER = 61.59453319166504 VAFACTOR= 0.0 ORIENTAT= 61.59453319166504 RA_APER = 149.5 DEC_APER= 69.8 EXPTIME = 1000.0 RDNOISE = 12.0 HISTORY Initialized Detector SCA01 with filter F158 HISTORY Adding items from catalogue catalog.txt HISTORY Adding 2 point sources HISTORY Added background of 0.7521 counts/s/detector pixel (0.7521 counts/s/over HISTORY sampled pixel) HISTORY Convolving SCA01 with PSF grid HISTORY Adding Poisson Noise with mean -0.007881569676101208 and standard deviat HISTORY ion 27.416147232055664 HISTORY Adding Read noise with mean -0.003449 and standard deviation 11.996563 HISTORY Adding Flatfield residual with mean 1.009788 and standard deviation 0.00 HISTORY 0677 HISTORY Adding Dark residual with mean 0.693514 and standard deviation 3.668765 HISTORY Adding Cosmic Ray residual with mean 0.279746 and standard deviation 3.4 HISTORY 87007

york-stsci commented 4 years ago

As such I can't reproduce this bug. For reference, my STIPS environment dictionary is as follows:

>>> import stips >>> stips.__env__dict__ { 'stips_version': '1.0.7.dev289', 'stips_data_location': '/Users/york/Data/Reference/stips_data-1.0.8', 'stips_data_version': '1.0.8', 'stips_grid_version': '1.0.7.dev', 'pandeia_version': '1.5', 'pandeia_data_location': '/Users/york/Data/Reference/pandeia_data-1.5_both', 'pandeia_data_version': '1.5', 'webbpsf_version': '0.9.0', 'webbpsf_data_location': '/Users/york/Data/Reference/webbpsf-data-0.9.0', 'webbpsf_data_version': '0.9.0', 'astropy_version': '4.0.1.post1', 'photutils_version': '0.7.2' }

I wonder if the issue is that there's an out-of-date python module that you're using?

benw1 commented 4 years ago

You did not set the pa to 0. I think that is necessary to reproduce the problem.

york-stsci commented 4 years ago

I now get:

XTENSION= 'IMAGE ' / Image extension BITPIX = -32 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 4096 NAXIS2 = 4096 PCOUNT = 0 / number of parameters GCOUNT = 1 / number of groups WCSAXES = 2 / Number of coordinate axes CRPIX1 = 2048.0 / Pixel coordinate of reference point CRPIX2 = 2048.0 / Pixel coordinate of reference point CDELT1 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point CDELT2 = 3.0555555555556E-05 / [deg] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection CRVAL1 = 149.5 / [deg] Coordinate value at reference point CRVAL2 = 69.8 / [deg] Coordinate value at reference point LONPOLE = 180.0 / [deg] Native longitude of celestial pole LATPOLE = 69.8 / [deg] Native latitude of celestial pole DATEREF = '1858-11-17' / ISO-8601 fiducial time MJDREFI = 0.0 / [d] MJD of fiducial time, integer part MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part RADESYS = 'ICRS' / Equatorial coordinate system EXTNAME = 'SCA01 ' / extension name DETECTOR= 'SCA01 ' FILTER = 'F158 ' DATE-OBS= '2020-09-30' TIME-OBS= 'Sep:09:1601490388' EQUINOX = 2000.0 PA_APER = 0.0 VAFACTOR= 0.0 ORIENTAT= 0.0 RA_APER = 149.5 DEC_APER= 69.8 EXPTIME = 1000.0 RDNOISE = 12.0 HISTORY Initialized Detector SCA01 with filter F158 HISTORY Adding items from catalogue catalog.txt HISTORY Adding 1 point sources HISTORY Added background of 0.7521 counts/s/detector pixel (0.7521 counts/s/over HISTORY sampled pixel) HISTORY Convolving SCA01 with PSF grid HISTORY Adding Poisson Noise with mean 0.007044952362775803 and standard deviati HISTORY on 27.426918029785156 HISTORY Adding Read noise with mean 0.003083 and standard deviation 12.001271 HISTORY Adding Flatfield residual with mean 1.009788 and standard deviation 0.00 HISTORY 0677 HISTORY Adding Dark residual with mean 0.693514 and standard deviation 3.668765 HISTORY Adding Cosmic Ray residual with mean 0.280163 and standard deviation 3.4 HISTORY 89845

So apparently astropy.wcs doesn't add in a PC or CD matrix if it would be trivial. If you want that to change, that's not something that STIPS can do. You'll need to file a bug against astropy.