Photometry

The following examples are adapted from Photometry.jl to show the same examples combined with AstroImages.jl. To learn how to measure background levels, perform aperture photometry, etc see the Photometry.jl documentation.

Background Estimation

From Photometry.jl:

Estimating backgrounds is an important step in performing photometry. Ideally, we could perfectly describe the background with a scalar value or with some distribution. Unfortunately, it's impossible for us to precisely separate the background and foreground signals. Here, we use mixture of robust statistical estimators and meshing to let us get the spatially varying background from an astronomical photo. Let's show an example Now let's try and estimate the background using estimate_background. First, we'll si gma-clip to try and remove the signals from the stars. Then, the background is broken down into boxes, in this case of size (50, 50). Within each box, the given statistical estimators get the background value and RMS. By default, we use SourceExtractorBackground and StdRMS. This creates a low-resolution image, which we then need to resize. We can accomplish this using an interpolator, by default a cubic-spline interpolator via ZoomInterpolator. The end result is a smooth estimate of the spatially varying background and background RMS.

using Photometry
using AstroImages
using Plots # optional, for implot functionality

# Download our image, courtesy of astropy
image = AstroImage(download("https://rawcdn.githack.com/astropy/photutils-datasets/8c97b4fa3a6c9e6ea072faeed2d49a20585658ba/data/M6707HH.fits"))

# sigma-clip
clipped = sigma_clip(image, 1, fill=NaN)

# get background and background rms with box-size (50, 50)
bkg, bkg_rms = estimate_background(clipped, 50)

imview(image)
imview(clipped)
imview(bkg)
imview(bkg_rms)

Or, if you have Plots loaded:

using Plots

AstroImages.set_clims!(Percent(99.5))
AstroImages.set_cmap!(:magma)
AstroImages.set_stretch!(identity)

plot(
    implot(image, title="Original"),
    implot(clipped, title="Sigma-Clipped"),
    implot(bkg, title="Background"),
    implot(bkg_rms, title="Background RMS"),
    layout=(2, 2)
)

We could apply a median filter, too, by specifying filter_size

# get background and background rms with box-size (50, 50) and filter_size (5, 5)
bkg_f, bkg_rms_f = estimate_background(clipped, 50, filter_size=5)

# plot
plot(
    implot(bkg, title="Unfiltered", ylabel="Background"),
    implot(bkg_f, title="Filtered"),
    implot(bkg_rms, ylabel="RMS"),
    implot(bkg_rms_f);
    layout=(2, 2),)

Now we can see our image after subtracting the filtered background and ready for Aperture Photometry!

subt = image .- bkg_f[axes(image)...]
clims = extrema(vcat(vec(image), vec(subt)))
plot(
    implot(image; title="Original", clims),
    implot(subt; title="Subtracted", clims),
    size=(1600,1000)
)

Source Extraction

From the background-subtracted image, we can detect all sources in the image:

# We specify the uncertainty in the pixel data. We'll set it equal to zero.
errs = zeros(axes(subt))
sources = extract_sources(PeakMesh(), subt, errs, true) # sort from brightest to darkest

Table with 3 columns and 60924 rows:
      x     y     value
    ┌────────────────────
 1  │ 255   226   9762.42
 2  │ 940   681   9661.5
 3  │ 219   924   9653.78
 4  │ 38    678   9647.93
 5  │ 245   85    9637.55
 6  │ 44    1001  9633.13
 7  │ 503   904   9627.34
 8  │ 819   775   9625.92
 9  │ 610   62    9617.79
 10 │ 133   110   9616.88
 11 │ 592   123   9614.86
 12 │ 1055  248   9612.46
 13 │ 424   86    9611.86
 14 │ 558   125   9585.93
 15 │ 125   1021  9582.11
 16 │ 112   637   9579.71
 17 │ 196   2     9579.49
 ⋮  │  ⋮     ⋮       ⋮

There's over 60,000 sources!

We'll define a circular apperture for each source:

aps = CircularAperture.(sources.x, sources.y, 6)[1:1000] # just brightest thousand point sources

1000-element Vector{CircularAperture{Int64}}:
 CircularAperture(255, 226, r=6)
 CircularAperture(940, 681, r=6)
 CircularAperture(219, 924, r=6)
 CircularAperture(38, 678, r=6)
 CircularAperture(245, 85, r=6)
 CircularAperture(44, 1001, r=6)
 CircularAperture(503, 904, r=6)
 CircularAperture(819, 775, r=6)
 CircularAperture(610, 62, r=6)
 CircularAperture(133, 110, r=6)
 ⋮
 CircularAperture(697, 975, r=6)
 CircularAperture(703, 599, r=6)
 CircularAperture(893, 878, r=6)
 CircularAperture(291, 457, r=6)
 CircularAperture(660, 855, r=6)
 CircularAperture(711, 920, r=6)
 CircularAperture(708, 271, r=6)
 CircularAperture(701, 273, r=6)
 CircularAperture(710, 273, r=6)

We can overplot them on our original image. The coordinate sytem used by the Photometry.jl plot recipes (but not the actual return values) doesn't match AstroImages, so we must transpose our image:

implot(subt', colorbar=false)
plot!(aps)

Measuring Photometry

Finally we can extract the source photometry

table = photometry(aps, subt)

Table with 3 columns and 1000 rows:
      xcenter  ycenter  aperture_sum
    ┌───────────────────────────────
 1  │ 255      226      20542.4
 2  │ 940      681      7125.3
 3  │ 219      924      14741.0
 4  │ 38       678      31056.4
 5  │ 245      85       -910.308
 6  │ 44       1001     2792.62
 7  │ 503      904      6954.07
 8  │ 819      775      -3136.61
 9  │ 610      62       4338.13
 10 │ 133      110      22219.0
 11 │ 592      123      18562.2
 12 │ 1055     248      3277.33
 13 │ 424      86       255299.0
 14 │ 558      125      6659.03
 15 │ 125      1021     11524.1
 16 │ 112      637      1349.78
 17 │ 196      2        -797.491
 ⋮  │    ⋮        ⋮          ⋮

And plot them:

scatter(
    table.xcenter,
    table.ycenter,
    aspectratio=1,
    marker_z=table.aperture_sum,
    markerstrokewidth=0,
    label="",
    framestyle=:box,
    background_inside=:black,
    color=:white
)