Commit bdca120b authored by Mohammad Akhlaghi's avatar Mohammad Akhlaghi
Browse files

MUSE and HST PSFs created, but images not matched

The MUSE PSFs were given as parameters to a Moffat function fit. So
they were created using Gnuastro's MakeProfiles. The HST PSF was also
generated from one of the bright stars in the UDF. It isn't the
brightest, but so far we haven't checked if it was saturated or
not. However, apparently, Gnuastro's convolve had problems in finding
the kernel to convolve with the HST image to get to the MUSE
resolution.

So for now, we have kept the rules (to work on them later), but they
are currently ignored in the pipeline (they aren't a prerequisite to
anything). Since the MUSE PSF is much wider than the HST, for now, we
will just convolve the HST image with the MUSE PSF. Later, when the
problem above is solved, we will find the proper kernel more
accurately.
parent a0b58e9f
......@@ -53,8 +53,8 @@ all: $(BSYM) file.txt
#
# Note that we cannot simply include `reproduce/src/*.mk'. Because the
# order of reading them into Make actually matters in some cases.
include $(foreach m, preparations download input-cutouts , \
reproduce/src/$(m).mk)
include $(foreach m, preparations download degrade-hst input-cutouts \
, reproduce/src/$(m).mk)
......
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input image:
hdu 0
hdu 0
hdu 0
hdu 0
hdu 0
hdu 0
hdu 0
# Output:
type float
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Operating mode:
spatial 0
frequency 1
makekernel 0
# Input:
hdu 0
kernel kernel.fits
khdu 0
# Output:
# Mesh grid:
meshsize 32
nch1 1
nch2 1
lastmeshfrac 0.6
fullconvolution 0
\ No newline at end of file
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input:
hdu 0
# Output:
\ No newline at end of file
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Operating mode:
......
# Reproduction pipeline
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input:
......
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input:
hdu 0
objhdu 1
clumphdu 2
skyhdu 3
stdhdu 4
zeropoint 0.0
skysubtracted 0
# Output:
nsigmag 1
intwidth 6
floatwidth 13
accuwidth 15
floatprecision 3
accuprecision 8
# Upper limit magnitude:
upnum 200
upsclipmultip 3
upsclipaccu 0.2
upnsigma 1
# Catalog columns:
# sn 1
# std 1
# sky 1
# rivernum 1
# riverave 1
# clumpsdec 1
# clumpsra 1
# clumpsy 1
# clumpsx 1
# clumpsmagnitude 1
# clumpsarea 1
# numclumps 1
# magnitude 1
# area 1
# dec 1
# ra 1
# y 1
# x 1
# idinhostobj 1
# hostobjid 1
id 1
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input:
hdu 0
# Output:
naxis1 1000
naxis2 1000
oversample 5
circumwidth 2
type float
# Profiles:
tunitinp 0
numrandom 10000
tolerance 0.01
zeropoint 0.00
prepforconv 0
xshift 0
yshift 0
# Catalog:
xcol 1
ycol 2
fcol 3
rcol 4
ncol 5
pcol 6
qcol 7
mcol 8
tcol 9
# WCS:
crpix1 1
crpix2 1
crval1 1
crval2 1
resolution 0.03
\ No newline at end of file
# Reproduction pipeline
nolog 1
onlydirconf 1
keepinputdir 1
#onlyversion 0.1
# Input:
hdu 0
khdu 0
skysubtracted 0
minbfrac 0.5
minnumfalse 100
# Output:
grownclumps 0
# Mesh grid:
smeshsize 32
lmeshsize 200
nch1 1
nch2 1
lastmeshfrac 0.51
numnearest 10
smoothwidth 3
fullconvolution 0
fullinterpolation 0
fullsmooth 0
# Detection:
mirrordist 1.5
minmodeq 0.49
qthresh 0.3
erode 2
erodengb 4
noerodequant 0.9331
opening 1
openingngb 8
sigclipmultip 3
sigcliptolerance 0.2
dthresh -0.1
detsnminarea 15
detsnhistnbins 0
detquant 0.95
dilate 3
# Segmentation
segsnminarea 25
keepmaxnearriver 0
segquant 0.95
clumpsnhistnbins 0
gthresh 0.5
minriverlength 15
objbordersn 1
\ No newline at end of file
# Configuration files for programs
mkprof = reproduce/config/gnuastro/astmkprof.conf
imgcrop = reproduce/config/gnuastro/astimgcrop.conf
imgwarp = reproduce/config/gnuastro/astimgwarp.conf
# The pixel scale of the input images in units of 0.01 arcseconds
#
# The reason it is in units of 0.01 arcseconds is that we want to use
# it in the filenames and this makes it easy.
hstpixelscale = 3
hst-star-ra = 53.132269
hst-star-dec = -27.782852
# The size of the PSFs from MUSE and HST to match.
psf-match-width = 21
# Only pixels above this multiple of the Sky standard deviation should
# be used in finding the center of the known star.
star-cat-threshold = 30
# Truncation radius for generating the MUSE PSF
#
# This is in units of the FWHM.
truncation-muse-psf = 3
# Broad-band photometry checks with MUSE generated broad-band images.
#
# Degrade the HST images to MUSE PSF and resolution
#
# Original author:
# Mohammad Akhlaghi <mohammad.akhlaghi@univ-lyon1.fr>
# Contributing author(s):
# Copyright (C) 2016, Mohammad Akhlaghi.
#
# This script is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This script is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# A copy of the GNU General Public License is available at
# <http://www.gnu.org/licenses/>.
# Top level directory for all cutouts.
dhstdir = $(BDIR)/degrade-hst
$(dhstdir): $(BDIR); mkdir $@
# Get the PSF information from the headers
all-muse-udf-psfs = $(foreach f, $(filters), $(MUSEINPUTS)/muse-udf-$(f).fits) \
$(foreach f, $(filters), $(MUSEINPUTS)/muse-udf10-$(f).fits)
$(dhstdir)/muse-psfs.txt: $(all-muse-udf-psfs) \
reproduce/config/internal/truncation-muse-psf.mk \
reproduce/config/internal/hst-pixel-scale.mk \
| $(dhstdir)
# Delete the output file because we will be appending to it,
# so if it already exists, there is going to be trouble.
rm -f $@
# Write the information for each PSF in the file. The
# MakeProfiles columns are defined in its configuration file.
for file in $^; do \
astheader $$file | grep FSF | \
awk 'BEGIN{print "#", "'$$file'"; bc=0; fc=0;} \
/FSF0.BET/ {beta[bc++]=$$NF} \
/FSF0.FWA/ {fwhm[fc++]=$$NF} \
END{ \
topixel=1/(0.01*$(hstpixelscale)); \
for(i=0;i<bc;++i) \
printf("1 0 0 1 %.3f %.3f 0 1 0 %.3f\n", \
fwhm[i]*topixel, beta[i], \
$(truncation-muse-psf)); \
}' >> $@; \
done
# Make the MUSE PSFs.
#
# The random number generator and seed must be fixed to get exactly
# the same results. For the seed, I just typed in some numbers to make
# a long integer. The important thing is that it be fixed, its
# absolute value is irrelevant.
muse-rawpsfs = $(dhstdir)/udf10-f850lp-rawpsf.fits
$(muse-rawpsfs): $(dhstdir)/muse-psfs.txt $(mkprof) | $(dhstdir)
# Generate all the PSFs using MakeProfiles. Note that
export GSL_RNG_TYPE="mt19937"; \
export GSL_RNG_SEED=129312398; \
astmkprof $< --individual --output=$(@D) --oversample=1 --nomerged
# Correct the names to a managable standard. You can see the
# commented lines in $(dhstdir)/muse-psfs.txt to check them.
mv $(@D)/0_muse-psfs.fits $(@D)/udf1-f606w-rawpsf.fits
mv $(@D)/1_muse-psfs.fits $(@D)/udf2-f606w-rawpsf.fits
mv $(@D)/2_muse-psfs.fits $(@D)/udf3-f606w-rawpsf.fits
mv $(@D)/3_muse-psfs.fits $(@D)/udf4-f606w-rawpsf.fits
mv $(@D)/4_muse-psfs.fits $(@D)/udf5-f606w-rawpsf.fits
mv $(@D)/5_muse-psfs.fits $(@D)/udf6-f606w-rawpsf.fits
mv $(@D)/6_muse-psfs.fits $(@D)/udf7-f606w-rawpsf.fits
mv $(@D)/7_muse-psfs.fits $(@D)/udf8-f606w-rawpsf.fits
mv $(@D)/8_muse-psfs.fits $(@D)/udf9-f606w-rawpsf.fits
mv $(@D)/9_muse-psfs.fits $(@D)/udf1-f775w-rawpsf.fits
mv $(@D)/10_muse-psfs.fits $(@D)/udf2-f775w-rawpsf.fits
mv $(@D)/11_muse-psfs.fits $(@D)/udf3-f775w-rawpsf.fits
mv $(@D)/12_muse-psfs.fits $(@D)/udf4-f775w-rawpsf.fits
mv $(@D)/13_muse-psfs.fits $(@D)/udf5-f775w-rawpsf.fits
mv $(@D)/14_muse-psfs.fits $(@D)/udf6-f775w-rawpsf.fits
mv $(@D)/15_muse-psfs.fits $(@D)/udf7-f775w-rawpsf.fits
mv $(@D)/16_muse-psfs.fits $(@D)/udf8-f775w-rawpsf.fits
mv $(@D)/17_muse-psfs.fits $(@D)/udf9-f775w-rawpsf.fits
mv $(@D)/18_muse-psfs.fits $(@D)/udf1-f814w-rawpsf.fits
mv $(@D)/19_muse-psfs.fits $(@D)/udf2-f814w-rawpsf.fits
mv $(@D)/20_muse-psfs.fits $(@D)/udf3-f814w-rawpsf.fits
mv $(@D)/21_muse-psfs.fits $(@D)/udf4-f814w-rawpsf.fits
mv $(@D)/22_muse-psfs.fits $(@D)/udf5-f814w-rawpsf.fits
mv $(@D)/23_muse-psfs.fits $(@D)/udf6-f814w-rawpsf.fits
mv $(@D)/24_muse-psfs.fits $(@D)/udf7-f814w-rawpsf.fits
mv $(@D)/25_muse-psfs.fits $(@D)/udf8-f814w-rawpsf.fits
mv $(@D)/26_muse-psfs.fits $(@D)/udf9-f814w-rawpsf.fits
mv $(@D)/27_muse-psfs.fits $(@D)/udf1-f850lp-rawpsf.fits
mv $(@D)/28_muse-psfs.fits $(@D)/udf2-f850lp-rawpsf.fits
mv $(@D)/29_muse-psfs.fits $(@D)/udf3-f850lp-rawpsf.fits
mv $(@D)/30_muse-psfs.fits $(@D)/udf4-f850lp-rawpsf.fits
mv $(@D)/31_muse-psfs.fits $(@D)/udf5-f850lp-rawpsf.fits
mv $(@D)/32_muse-psfs.fits $(@D)/udf6-f850lp-rawpsf.fits
mv $(@D)/33_muse-psfs.fits $(@D)/udf7-f850lp-rawpsf.fits
mv $(@D)/34_muse-psfs.fits $(@D)/udf8-f850lp-rawpsf.fits
mv $(@D)/35_muse-psfs.fits $(@D)/udf9-f850lp-rawpsf.fits
mv $(@D)/36_muse-psfs.fits $(@D)/udf10-f606w-rawpsf.fits
mv $(@D)/37_muse-psfs.fits $(@D)/udf10-f775w-rawpsf.fits
mv $(@D)/38_muse-psfs.fits $(@D)/udf10-f814w-rawpsf.fits
mv $(@D)/39_muse-psfs.fits $(@D)/udf10-f850lp-rawpsf.fits
# Cutout the MUSE PSFs to the proper size
#
# The size of the images created by MakeProfiles was based on the PSF
# FWHM. But to match with the HST PSF, we need a different size (set
# in the reproduction pipeline configurations). So here, we will crop
# the desired region from the raw psfs.
muse-psfs = $(foreach uid, 1 2 3 4 5 6 7 8 9 10, \
$(foreach f, $(filters), $(dhstdir)/udf$(uid)-$(f)-psf.fits) )
$(muse-psfs): %-psf.fits: $(muse-rawpsfs) \
reproduce/config/internal/psf-match-width.mk
# First do a sanity check, since the size is fixed. Note that
# the image sizes are always an odd number, and since AWK does
# all internal arithmetic in floating point, we will just
# subtract one after dividing by two to get the central point.
xyc=$$(astheader $*-rawpsf.fits | grep NAXIS \
| awk 'NR==2{xw=$$3} NR==3{yw=$$3} \
END{printf "--xc=%d --yc=%d", xw/2+1, yw/2+1}'); \
astimgcrop $*-rawpsf.fits $$xyc -o$@ --zeroisnotblank \
--iwidth=$(psf-match-width)
# HST PSF
#
# Crop out one star from the HST image, center it and then use it for
# the HST PSF.
xdf-prefix = $(XDF)/hlsp_xdf_hst_acswfc-$(hstpixelscale)0mas_hudf_
hst-star-cats = $(foreach f, $(filters), $(dhstdir)/star-$(f).txt)
$(hst-star-cats): $(dhstdir)/star-%.txt: $(xdf-prefix)%_v1_sci.fits \
reproduce/config/internal/hst-star-position.mk \
reproduce/config/internal/star-cat-threshold.mk | $(dhstdir)
# Crop the region in the vicinity of the star in a large
# enough box that the Sky can be accurate calculated. A width
# of 15 arcseconds corresponds to a 500x500 box.
astimgcrop --ra=$(hst-star-ra) --dec=$(hst-star-dec) --wwidth=15 \
$< -o$(@D)/$*-l.fits
# Run NoiseChisel on the image
astnoisechisel $(@D)/$*-l.fits
# Make a catalog. Note that since we want the center of the
# star to be measured accurately, we use a very large
# threshold. We are just getting the area
astmkcatalog $(@D)/$*-l_labeled.fits \
--area --y --x --threshold=$(star-cat-threshold)
# The star should be in the central position: the object
# closest to the coordinates (250,250). Since this is a bright
# star, we don't expect any star within a few pixels.
awk '!/^#/ && $$2>248 && $$2<252 && $$3>248 && $$3<252 \
{print $$2, $$3}' $(@D)/$*-l_labeled_c.txt > $@
# Clean up and set final name
rm $(@D)/$*-l_labeled.fits
rm $(@D)/$*-l_labeled_o.txt $(@D)/$*-l_labeled_c.txt
# Using the central position of the star calculated from the catalog,
# warp the image into a grid with the star center exactly in the
# center of a pixel
hst-centered-star = $(foreach f, $(filters), $(dhstdir)/star-$(f).fits)
$(hst-centered-star): $(dhstdir)/star-%.fits: $(dhstdir)/star-%.txt \
reproduce/config/internal/psf-match-width.mk
# Generate the proper options for ImageWarp, then run it. As
# an example, when the position along an axis is 250.479, we
# need to translate the grid by 250-250.479 = -0.479 so the
# center of the star moves to the center of the pixel.
iwopts=$$(awk '{printf("--translate=%.3f,%.3f", 250-$$1, 250-$$2)}' \
$<); \
astimgwarp $(@D)/$*-l.fits $$iwopts -o$(@D)/$*-w.fits
# Crop the centeral few pixels of the image to exactly
# pinpoint the center.
astimgcrop --xc=250 --yc=250 $(@D)/$*-w.fits --iwidth=5
astconvertt $(@D)/$*-w_crop.fits -o$(@D)/$*-c.txt
rm $(@D)/$*-w_crop.fits
# Identify the peak pixel.
xyc=$$(awk 'BEGIN{max=-9e15} \
!/^#/{++rc; \
for(cc=1;cc<=NF;++cc) \
if($$cc>max) {max=$$cc; xc=cc; yc=rc;} \
} \
END{printf "--xc=%d --yc=%d", 250+xc-3, 250+yc-3}' \
$(@D)/$*-c.txt); \
astimgcrop $$xyc $(@D)/$*-w.fits --iwidth=$(psf-match-width) \
-o$(@D)/$*-cc.txt
# Clean up
rm $(@D)/$*-c.txt $(@D)/$*-w.fits
# Kernels to match HST and MUSE
#
# Using the PSF images for MUSE and HST, now we want to find the
# kernel that the HST images should be convolved with to have similar
# PSF to MUSE.
kernels = $(foreach uid, 1 2 3 4 5 6 7 8 9 10, \
$(foreach f, $(filters), $(dhstdir)/udf$(uid)-$(f)-ker.fits) )
$(kernels): $(dhstdir)/udf%-ker.fits: $(dhstdir)/udf%-psf.fits
# Get the HST image name and find the kernel
filter=$(word 2, $(subst -, , $*)); \
astconvolve --kernel=$< $(dhstdir)/star-$$filter".fits" --makekernel=8 \
--viewfreqsteps
file.txt: $(muse-psfs)
echo; echo "-- REACHED THE END ---"; echo
......@@ -113,10 +113,3 @@ udf-muse-cutouts = $(udf1-muse-cutouts) $(udf2-muse-cutouts) \
$(udf5-muse-cutouts) $(udf6-muse-cutouts) \
$(udf7-muse-cutouts) $(udf8-muse-cutouts) \
$(udf9-muse-cutouts) $(udf10-muse-cutouts)
file.txt: $(udf-muse-cutouts)
echo
echo "end of input-cutouts"
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment