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(
    (
        band=:L,
        image=AstroImages.recenter(image), platescale=9.971,
        epoch=mjd("2021")
    ),
)
OctofitterImages.ImageLikelihood Table with 6 columns and 1 row:
     band  image                 platescale  epoch    contrast              ⋯
   ┌─────────────────────────────────────────────────────────────────────────
 1 │ L     [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 Visual{FixedPosition} begin
    sep ~ Uniform(0, 2000)
    pa ~ Uniform(0,2pi)
    # Contrast ratio
    L ~ Uniform(0, 1)
end imglike

@system sys begin
    plx = 24.4620
end b

model = Octofitter.LogDensityModel(sys, verbosity=4)
LogDensityModel for System sys of dimension 3 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.

model.starting_points = model.link.([
    [1704, deg2rad(70.63), 1e-4]
])
chain = octofit(model, iterations=10000)
Chains MCMC chain (10000×17×1 Array{Float64, 3}):

Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 4.82 seconds
Compute duration  = 4.82 seconds
parameters        = plx, b_sep, b_pa, b_L
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      rhat ⋯
      Symbol     Float64   Float64   Float64     Float64     Float64   Float64 ⋯

         plx     24.4620    0.0000    0.0000         NaN         NaN       NaN ⋯
       b_sep   1706.7958    7.5897    1.2643     71.7172    100.9916    1.0004 ⋯
        b_pa      1.2331    0.0020    0.0000   5411.2530   6804.2113    1.0000 ⋯
         b_L      0.0001    0.0000    0.0000   2892.0153   3164.0368    0.9999 ⋯
                                                                1 column omitted

Quantiles
  parameters        2.5%       25.0%       50.0%       75.0%       97.5%
      Symbol     Float64     Float64     Float64     Float64     Float64

         plx     24.4620     24.4620     24.4620     24.4620     24.4620
       b_sep   1698.7710   1702.0877   1704.1023   1707.9655   1725.5502
        b_pa      1.2296      1.2315      1.2329      1.2345      1.2370
         b_L      0.0001      0.0001      0.0001      0.0001      0.0001

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
model.starting_points = nothing # reset starting points
chain, pt = octofit_pigeons(model, n_rounds=11)
(chain = MCMC chain (2048×8×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))
The median separation is 1703.8388560973906

Visualize

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