fluximage & friends w/ random binsize

[kglotfelty]

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]

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]

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]

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]

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!