fluximage & friends w/ random binsize

kglotfelty commented 1 year ago

Ref: cxc helpdesk ticket 24414.

There is a precision issue when a users uses a binsize with 100's digits and beyond.

% download_chandra_obsid 4425 --quiet
% cd 4425
% fluximage . out=boo bin=2.35 clob+ psfecf=0.5 mode=h clob+
Running fluximage
Version: 04 November 2021

Found ./primary/acisf04425N005_evt2.fits.gz
Using event file ./primary/acisf04425N005_evt2.fits.gz
Using CSC ACIS broad science energy band.
Aspect solution primary/pcadf04425_000N001_asol1.fits found.
Bad-pixel file primary/acisf04425_000N005_bpix1.fits.gz found.
Mask file secondary/acisf04425_000N005_msk1.fits found.

The output images will have 406 by 326 pixels, pixel size of 1.1562000000000001 arcsec,
    and cover x=3405.65:4359.75:2.35,y=3520.8:4286.900000000001:2.35.

Running tasks in parallel with 4 processors.
Creating aspect histogram for obsid 4425
Creating instrument map for obsid 4425
Creating exposure map for obsid 4425
Thresholding data for obsid 4425
Exposure-correcting image for obsid 4425
Creating PSF map for obsid 4425
# mkpsfmap (CIAO 4.15): ERROR: Cannot load infile 'boo_broad_thresh.img[sky=MASK(boo_broad_thresh.expmap)]'

# mkpsfmap (CIAO 4.15): ERROR: Cannot open input file 'boo_broad_thresh.img[sky=MASK(boo_broad_thresh.expmap)]'

Exception ignored in: <function _TemporaryFileCloser.__del__ at 0x7fcaeec39310>
Traceback (most recent call last):
  File "/export/ciao_from_source/ciao-4.15/ots/lib/python3.9/tempfile.py", line 445, in __del__
    self.close()
  File "/export/ciao_from_source/ciao-4.15/ots/lib/python3.9/tempfile.py", line 441, in close
    unlink(self.name)
FileNotFoundError: [Errno 2] No such file or directory: '/tmp/tmpy_7wj_em.psfmap'
# fluximage (04 November 2021): ERROR Unable to run mkpsfmap with arguments: infile=boo_broad_thresh.img[sky=MASK(boo_broad_thresh.expmap)] outfile=/tmp/tmpy_7wj_em.psfmap[PSFMAP] energy=2.3 spectrum= ecf=0.5 mode=h clobber=yes

The problem is that the WCS in the exposure map and the image no longer match so trying to use the expmap as a mask fails:

% dmlist 'boo_broad_thresh.img[sky=MASK(boo_broad_thresh.expmap)]' blocks

Failed to open virtual file boo_broad_thresh.img[sky=MASK(boo_broad_thresh.expmap)]
# dmlist (CIAO 4.15): Mask coordinate failure, mismatch 'LINEAR' offset for 'SKY'.

% dmlist boo_broad_thresh.expmap cols | egrep -i '\(X\)|\(Y\)'
   1   1,2    SKY(X) = (+3406.8750) +(+2.350 )* ((#1)-(+1.0))
                 (Y)   (+3521.9749)  (+2.3498)  ((#2) (+1.0))
   1   1,2    EQPOS(RA ) = (+173.2652)[deg] +TAN[(-0.000136667)* (SKY(X)-(+4096.50))]
                   (DEC)   (+25.9069 )           (+0.000136667)  (   (Y) (+4096.50)) 
% dmlist boo_broad_thresh.img cols | egrep -i '\(X\)|\(Y\)'
   1   1,2    sky(x) = (+3405.650) +(+2.350)* ((#1)-(+0.50))
                 (y)   (+3520.80 )  (+2.350)  ((#2) (+0.50))
   1   1,2    EQPOS(RA ) = (+173.2652)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+25.9069 )           (+0.000136667)  (   (y) (+4096.50))

we can see that the y-axis pixel size in the exposure map is diff than the x-axis, and diff from the image.

DougBurke commented 1 year ago

So, why do the WCS's no-longer match? Why is the exposure map using a Y size of 2.3498????? What foul and wicky deed did I make the code do this time?

DougBurke commented 1 year ago

So, we get, from a run with verbose=2

Running tool skyfov using:
>>> skyfov primary/acisf04425N005_evt2.fits.gz /tmp/tmpcnwrhsyg.fov mskfile=secondary/acisf04425_000N005_msk1.fits aspect="primary/pcadf04425_000N001_asol1.fits[@primary/acisf04425N005_evt2.fits.gz]" method=convexhull clobber=yes
The output images will have 406 by 326 pixels, pixel size of 1.1562000000000001 arcsec,
    and cover x=3405.65:4359.75:2.35,y=3520.8:4286.900000000001:2.35.

This y upper limit could be an issue. It is used - e.g. in

Running tasks in parallel with 4 processors.
Creating aspect histogram for obsid 4425
>> Running asphist
>>   args: ['infile=primary/pcadf04425_000N001_asol1.fits', 'outfile=boo_7.asphist', 'evtfile=primary/acisf04425N005_evt2.fits.gz[ccd_id=7]', 'dtffile=', 'max_bin=10000', 'res_xy=0.5781000000000001', 'mode=h', 'clobber=yes', 'verbose=1']
>> Running dmcopy
>>   args: ['infile=primary/acisf04425N005_evt2.fits.gz[bin x=3405.65:4359.75:2.35,y=3520.8:4286.900000000001:2.35][energy=500.0:7000.0][opt type=i4]', 'outfile=/tmp/tmp5f172v4g.img', 'mode=h', 'clobber=yes', 'verbose=1']

Later on we end up with

Creating exposure map for obsid 4425
>> Running get_sky_limits
>>   args: ['image=boo_broad.img', 'verbose=0', 'mode=h']
>> Running mkexpmap
>>   args: ['outfile=boo_7_broad.expmap', 'instmapfile=boo_7_broad.instmap', 'asphistfile=boo_7.asphist[asphist]', 'xygrid=3405.7:4359.8:#406,3520.8:4289.2:#327', 'normalize=no', 'useavgaspect=no', 'mode=h', 'clobber=yes', 'verbose=1']

Note here I've used "#xxx" to define the number of bins. We can compare it to

bin x=3405.65:4359.75:2.35,y=3520.8:4286.900000000001:2.35

so it is slightly different.

% dmlist boo_broad_thresh.img subspace | grep -i sky
   1 sky                  [ 1] x                   3405.650:     4359.750 
   1 sky                  [ 2] y                   3520.80:     4289.250 
% dmlist boo_broad_thresh.expmap subspace | grep -i sky
   1 SKY                  [ 1] X                   3405.70:     4359.80 
   1 SKY                  [ 2] Y                   3520.80:     4289.20

and

% dmlist boo_broad_thresh.img cols

--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------

ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[406,327]              Int4(406x327)  -                    

--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------

Group# Axis# 
   1   1,2    sky(x) = (+3405.650) +(+2.350)* ((#1)-(+0.50))
                 (y)   (+3520.80 )  (+2.350)  ((#2) (+0.50))

--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------

Group# Axis# 
   1   1,2    EQPOS(RA ) = (+173.2652)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+25.9069 )           (+0.000136667)  (   (y) (+4096.50))

and

% dmlist boo_broad_thresh.expmap cols

--------------------------------------------------------------------------------
Columns for Image Block EXPMAP
--------------------------------------------------------------------------------

ColNo  Name                 Unit        Type             Range
   1   EXPMAP[406,327]      cm**2 s      Real4(406x327) -Inf:+Inf            

--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EXPMAP
--------------------------------------------------------------------------------

Group# Axis# 
   1   1,2    SKY(X) = (+3406.8750) +(+2.350 )* ((#1)-(+1.0))
                 (Y)   (+3521.9749)  (+2.3498)  ((#2) (+1.0))

--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EXPMAP
--------------------------------------------------------------------------------

Group# Axis# 
   1   1,2    EQPOS(RA ) = (+173.2652)[deg] +TAN[(-0.000136667)* (SKY(X)-(+4096.50))]
                   (DEC)   (+25.9069 )           (+0.000136667)  (   (Y) (+4096.50))

Oh joy.

kglotfelty commented 1 year ago

I think the gist is that when the grids you use for the exposure maps that are actually used are limited to one digit of precision so

% dmhistory boo_broad_thresh.expmap mkexpmap
mkexpmap asphistfile="boo_7.asphist[asphist]" outfile="boo_7_broad.expmap" 
instmapfile="boo_7_broad.instmap" 
xygrid="3405.7:4359.8:#406,3520.8:4289.2:#327" 
normalize="no" useavgaspect="no" geompar="geom" verbose="0" clobber="yes"

I think when I also tried to use the xygrid parameter that it still translates it to use the #pixels format rather than using the binsize which surprised me.

DougBurke commented 1 year ago

I think I used #foo because I'd seem problems at other times when not doing it. Thoughts include

punting, and just manually re-writing the WCS of the exposure map - but this is hard to actually do
dropping the accuracy limit, but I'm worried what other consequences there are
try to pass around the "correct" xygrid (which is probably just the previous issue, but maybe with it's own fun)
add some validation so we can error out when we know there's a difference (not at all useful to the user but would be internaly to try and identify when this happens).

This all strikes me as something that can't be addressed quickly!

cxcsds / ciao-contrib

fluximage & friends w/ random binsize #710