Re: Using astscript-psf-select-stars on HST flc images

[Mohammad Akhlaghi]

Hi again Alex,

I was able to track down the problem: indeed, there are 4 WCSs within 
the header file! You can see with the second command below, where each 
set of WCS parameters starts with a 'WCSNAME' (I am including the 'wget' 
for future reference):

wget -Oj9en0pnaq_flc.fits 
'https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/j9en0pnaq_flc.fits' 


$ astfits j9en0pnaq_flc.fits -h1 | grep WCSNAME
WCSNAMEO= 'OPUS    '           / Coordinate system title 

WCSNAME = 'IDC_4bb1536oj-FIT_REL_GAIAeDR3' / Coordinate system title 

WCSNAMEA= 'IDC_4bb1536oj'      / Coordinate system title 

WCSNAMEB= 'IDC_4bb1536oj-GSC240-1' / Coordinate system title

I'll have to consult Mark (author of WCSLIB) on this, to see if it is a 
non-standard issue with the keywords, a bug in WCSLIB or maybe a problem 
in the way we use WCSLIB in Gnuastro.

But until then, here is the solution: use the Crop program's '--hendwcs' 
option! With this option, Crop will only use the keywords before the 
given value. I defined it specifically for such issues I have confronted 
with WCS here and there! You can get the counter of each keyword with 
this command:

$ astfits j9en0pnaq_flc.fits -h1 | cat -n

You can see that the main 'WCSNAME' is on line 185. So if you run Crop 
like the following non of the pixels will be touched (because we are 
setting '--section=:,:' to use all the pixels), but the header keywords 
will be clean:

$ astcrop j9en0pnaq_flc.fits --mode=img --section=:,:  --hendwcs=186 
-oj9en0pnaq-h1.fits

The output of this matches the raw image nicely when I fit by WCS in DS9 
perfectly (when I set "tile" to "single" and "blink"). But now, no 
Gnuastro program will complain about the WCS and their outputs will have 
a WCS :-).

Just some extra tips:

TIP1: On many raw exposures, the bias/dark isn't estimated/subtracted 
nicely. So the channels are usually still visible. To deal with this, 
NoiseChisel can do its interpolations (for the thresholds, and sky) on 
different channels separately (while the detection map will cover the 
whole image). In this particular case, since each HDU has two channels 
(aligned along the horizontal), you can run NoiseChisel like the following:

$ astnoisechisel j9en0pnaq-h1.fits --numchannels=2,1 -oj9en0pnaq-h1-det.fits

When you later look at the 'SKY' extension in DS9, you will clearly see 
different values in each channel. You can quantitatively see this 
differing distribution in the ASCII histogram of Segment on the 'SKY' 
extension also (where there are two distributions that are clearly 
separate):

$ aststatistics j9en0pnaq-h1-det.fits -hSKY
Statistics (GNU Astronomy Utilities) 0.19
-------
Input: j9en0pnaq-h1-det.fits (hdu: SKY)
Unit: ELECTRONS
-------
  Number of elements:                      8388608
  Minimum:                                 3.403336e+01
  Maximum:                                 3.592985e+01
  Median:                                  3.500706e+01
  Mean:                                    34.98684029
  Standard deviation:                      0.7641997643
-------
Histogram:
 |                                                                 *
 |                                                                 *
 |                                                                 *
 |                                                                 *
 |         *                                                       *
 |         *                                                       *
 |*        *                                                       *
 |*        *                                                   * * *
 |*      ****                                             *    *****
 |* ************                                          **  ******
 |******************                                    ****************
 |----------------------------------------------------------------------



TIP2: as far as I remember, these raw images haven't yet been corrected 
for the tilting of the HST detector compared to the focal plane. This 
tilting will affect the extended PSF and it is best to correct that 
tilting first, then estimate the extended PSF of each exposure. You can 
just mask the cosmic rays and subtract the sky before warping.

I am curious how the raw WCS parameters work with the tilting. Can you 
compare the output of Gnuastro's Warp with the default HST pipeline's 
warping tool (I think it was DrizzlePac, see TIP3 below on the 
'-nocs.fits' image):

astwarp j9en0pnaq-h1-det-nocs.fits --output=aligned.

Gnuastro's Warp is designed to preserve the diffuse signal by using 
area-based resampling, and we have a full tutorial on its optimal usage:

https://www.gnu.org/software/gnuastro/manual/html_node/Moire-pattern-and-its-correction.html

TIP3: For masking cosmic rays, you can try this Arithmetic command to do 
the job pretty fast. On first look, its results look good, but you can 
fine-tune the constants if more closer inspection shows incompleteness 
or impurity in masking/detecting the cosmic rays. It is a good exercise 
in reverse polish notation ;-). If it isn't clear, let us know in 
another post.

hdu=1
in=j9en0pnaq-h1-det.fits
out=j9en0pnaq-h1-det-nocs.fits
astarithmetic $in -h$hdu set-i i i 0 eq nan where set-i \
              5 5 i filter-median set-f i i f - f / 5 gt \
              2 dilate nan where --output=$out

Cheers,
Mohammad