GitHub

Extracting Traditional Photometry and Astrometry

Though not its primary purpose, you can use Octofitter to extract traditional astrometry and photometry from one or more images. This uses the functionality in the Fit Orbits to Images tutorial, but with a much simpler model.

Instead of fitting an entire orbit, we will simply fit an X / Y position and brightness.

Start by loading your images:

using Octofitter
using OctofitterImages
using Distributions
using Pigeons
using AstroImages
using CairoMakie

# Load individual iamges
# image1 = load("image1.fits")
# image2 = load("image2.fits")

# Or slices from a cube:
# cube = load("cube1.fits")
# image1 = cube[:,:,1]

# Download sample images from GitHub
download(
    "https://zenodo.org/records/6823071/files/HR8799.2021.fits?download=1",
    "HR8799-2021.fits"
)

# Or multi-extension FITS (this example)
image = AstroImages.load("HR8799-2021.fits")
Example block output

You can preview the image using imview from AstroImages:

imview(image)
Example block output

Note that to accurately extract astrometry and photometry, the input image should have already been convolved with the star or planet point spread function. If this isn't available, a convolution by a Gaussian or Airy disk might be an acceptable approximation.

Build the model

First, we create a table of our image data that will be attached to the Planet:

imglike = ImageLikelihood(
    Table(
        image=[AstroImages.recenter(image)],
        platescale=[9.971],
        epoch=[mjd("2021")]
    ),
    variables=@variables begin
        # Planet flux in image units -- could be contrast, mags, Jy, or arb. as long as it's consistent with the units of the data you provide
        flux ~ Uniform(0, 1)
        # The following are optional parameters for marginalizing over instrument systematics:
        # Platescale uncertainty multiplier [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
        platescale = 1.0
        # North angle offset in radians [could use: northangle ~ Normal(0, deg2rad(1))]
        northangle = 0.0
    end
)
OctofitterImages.ImageLikelihood Table with 5 columns and 1 row:
     image                 platescale  epoch    contrast              ⋯
   ┌───────────────────────────────────────────────────────────────────
 1 │ [NaN NaN NaN NaN Na…  9.971       59215.0  218-element extrapo…  ⋯

Note that you can also supply a contrast curve or map directly. If not provided, a simple contrast curve will be calculated directly from the data.

Next create the simplest possible model of 2D position, plus a contrast variable matching the band name used in the ImageLikelihood above:

planet_b = Planet(
    name="b",
    basis=Visual{Octofitter.FixedPosition},
    likelihoods=[imglike],
    variables=@variables begin
        sep ~ Uniform(0, 2000)
        pa ~ Uniform(0,2pi)
        # Contrast ratio
    end
)

sys = System(
    name="sys",
    companions=[planet_b],
    likelihoods=[],
    variables=@variables begin
        plx = 24.4620
    end
)

model = Octofitter.LogDensityModel(sys, verbosity=4)
LogDensityModel for System sys of dimension 3 and 1 epochs with fields .ℓπcallback and .∇ℓπcallback

Sample from the model (locally)

If you already know where the planet is and you only want to extract astrometry from that known location, you can specify a starting point and use hamiltonian monte carlo as follows. This will be very very fast.

initialize!(model, (;
    planets=(;
        b=(;
            sep=1704,
            pa=deg2rad(70.63),
        )
    )
))
chain = octofit(model, iterations=10000)
Chains MCMC chain (10000×19×1 Array{Float64, 3}):

Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 3.56 seconds
Compute duration  = 3.56 seconds
parameters        = plx, b_sep, b_pa, b_images_flux, b_images_platescale, b_images_northangle
internals         = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error

Summary Statistics
           parameters        mean       std      mcse    ess_bulk    ess_tail  ⋯
               Symbol     Float64   Float64   Float64     Float64     Float64  ⋯

                  plx     24.4620    0.0000    0.0000         NaN         NaN  ⋯
                b_sep   1703.8117    3.1609    0.0655   2647.7935   3021.2836  ⋯
                 b_pa      1.2330    0.0021    0.0000   4363.7940   5808.0516  ⋯
        b_images_flux      0.0001    0.0000    0.0000   3136.7033   3492.3443  ⋯
  b_images_platescale      1.0000    0.0000       NaN         NaN         NaN  ⋯
  b_images_northangle      0.0000    0.0000       NaN         NaN         NaN  ⋯
                                                               2 columns omitted

Quantiles
           parameters        2.5%       25.0%       50.0%       75.0%       97 ⋯
               Symbol     Float64     Float64     Float64     Float64     Floa ⋯

                  plx     24.4620     24.4620     24.4620     24.4620     24.4 ⋯
                b_sep   1698.6126   1701.6947   1703.3001   1705.6004   1710.9 ⋯
                 b_pa      1.2297      1.2314      1.2327      1.2345      1.2 ⋯
        b_images_flux      0.0001      0.0001      0.0001      0.0001      0.0 ⋯
  b_images_platescale      1.0000      1.0000      1.0000      1.0000      1.0 ⋯
  b_images_northangle      0.0000      0.0000      0.0000      0.0000      0.0 ⋯
                                                                1 column omitted

Sample from the model (globally)

You could also try sampling across the entire image, without necessarily specifying a starting position. Note that if there are multiple candidates, taking the naive mean and standard deviation will average across all planets.

using Pigeons
initialize!(model)
chain, pt = octofit_pigeons(model, n_rounds=11)
(chain = MCMC chain (2048×10×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Access results

samples_sep = chain[:b_sep]
samples_pa = chain[:b_pa]
println("The median separation is ", median(samples_sep))

flux = chain[:b_images_flux]
println("The flux is ", mean(flux), " ± ", std(flux))
println("The \"SNR\" is ", mean(flux)/std(flux))
The median separation is 1703.822418142548
The flux is 7.744552380090894e-5 ± 6.420402598667364e-6
The "SNR" is 12.062409266512933

Visualize

using CairoMakie, PairPlots
octocorner(model,chain)
Example block output