Examples
Plotting
We have recipes for all our aperture types, so you can easily create overlays on your images.
using Photometry
using Plots
plot(CircularAperture(2, 3, 4), c=1, xlims=(-1, 12), ylims=(0, 9))
plot!(CircularAnnulus(5, 5, 2.1, 3), c=2)
plot!(EllipticalAperture(0, 0, 10, 1, 32), c=3)
plot!(EllipticalAnnulus(5, 5, 4, 5, 2, -32), c=4)
plot!(RectangularAperture(0, 0, 4, 4, 4), c=5)
plot!(RectangularAnnulus(5, 1, 3, 4, 4, 4), c=6)Simple Stars
Here is an example where we will find aperture fluxes for stars from M67. The dataset is provided as part of the astropy/photutils-datasets repository.
Let's start by downloading and showing our image
using Photometry
using Plots
using FITSIO
# Load data in
url = "https://rawcdn.githack.com/astropy/photutils-datasets/8c97b4fa3a6c9e6ea072faeed2d49a20585658ba/data/M6707HH.fits"
hdu = FITS(download(url))
chunk = read(hdu[1], 81:155, 71:150)
# Plot
function imshow(image; kwargs...)
xs, ys = axes(image)
data = transpose(image)
heatmap(xs, ys, data; aspect_ratio=1, xlim=extrema(xs), ylim=extrema(ys), kwargs...)
end
imshow(chunk)Now let's add some apertures!
positions = [
[47.5 , 67.5],
[29.5 , 62.5],
[23.5 , 48.5],
[17.5 , 29.5],
[13.25, 10.5],
[65.5 , 14.0]
]
radii = [3, 3, 2.7, 2, 2.7, 3]
aps = CircularAperture.(positions, radii)6-element Vector{CircularAperture{Float64}}:
CircularAperture(47.5, 67.5, r=3.0)
CircularAperture(29.5, 62.5, r=3.0)
CircularAperture(23.5, 48.5, r=2.7)
CircularAperture(17.5, 29.5, r=2.0)
CircularAperture(13.25, 10.5, r=2.7)
CircularAperture(65.5, 14.0, r=3.0)now let's plot them up
imshow(chunk)
plot!(aps, c=:white)and finally let's get our output table for the photometry
table = photometry(aps, chunk)Table with 3 columns and 6 rows:
xcenter ycenter aperture_sum
┌───────────────────────────────
1 │ 47.5 67.5 2.48267e5
2 │ 29.5 62.5 2.25989e5
3 │ 23.5 48.5 1.49979e5
4 │ 17.5 29.5 72189.4
5 │ 13.25 10.5 1.48118e5
6 │ 65.5 14.0 2.02803e5Stars with Spatial Background Subtraction and PSF Fitting
This example will be the same as Simple Stars but will add background estimation with BackgroundMeshes.jl and PSF fitting with PSFModels.jl.
# `sigma_clip` and `estimate_background` are reexported from BackgroundMeshes.jl for convenience
clipped = sigma_clip(chunk, 1, fill=NaN)
# Estimate 2D spatial background using boxes of size (5, 5)
bkg, bkg_rms = estimate_background(clipped, 5)
plot(
imshow(chunk, title="Original"),
imshow(clipped, title="Sigma-Clipped"),
imshow(bkg, title="Background"),
imshow(bkg_rms, title="Background RMS");
layout=(2, 2), size=(600, 600), ticks=false
)Now, using the same apertures, let's find the output using the background-subtracted image
plot(
imshow(chunk, title="Original"),
imshow(chunk .- bkg, title="Subtracted");
layout=2, size=(600, 260), ticks=false, colorbar=false
)
plot!(aps, c=:white, subplot=1)
plot!(aps, c=:white, subplot=2)using PSFModels
function fit_psf(img_ap)
# Normalize
psf_data = collect(Float32, img_ap)
psf_data ./= maximum(psf_data)
# Set params
y, x = Tuple(argmax(psf_data))
fwhm = 5.0
params = (; x, y, fwhm)
# Fit
psf_params, psf_model = PSFModels.fit(gaussian, params, psf_data; x_abstol = 2e-6)
# Could also return a Tuple to display more information.
# Just returning fitted FWHM here for simplicity.
return psf_params.fwhm
end
table = photometry(aps, chunk .- bkg, bkg_rms; f = fit_psf)Table with 5 columns and 6 rows:
xcenter ycenter aperture_sum aperture_sum_err aperture_f
┌─────────────────────────────────────────────────────────────
1 │ 47.5 67.5 2.13534e5 431.48 4.19789
2 │ 29.5 62.5 114217.0 887.37 3.46587
3 │ 23.5 48.5 59230.7 1061.28 2.58725
4 │ 17.5 29.5 23159.6 697.556 2.32789
5 │ 13.25 10.5 54638.1 1048.67 2.51242
6 │ 65.5 14.0 91179.1 1168.71 3.09683