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.52 seconds
Compute duration  = 3.52 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   691.0883   12.3236    0.4215    963.1591   1701.2175   ⋯
                 b_pa     4.1085    0.0060    0.0001   3821.8769   4608.7875   ⋯
        b_images_flux     0.0002    0.0000    0.0000   2101.0497   1946.1877   ⋯
  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.5%
               Symbol    Float64    Float64    Float64    Float64    Float64

                  plx    24.4620    24.4620    24.4620    24.4620    24.4620
                b_sep   666.6333   684.6601   689.5792   698.9191   715.1637
                 b_pa     4.0964     4.1045     4.1088     4.1124     4.1200
        b_images_flux     0.0001     0.0002     0.0002     0.0002     0.0002
  b_images_platescale     1.0000     1.0000     1.0000     1.0000     1.0000
  b_images_northangle     0.0000     0.0000     0.0000     0.0000     0.0000

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.8169793919144
The flux is 7.768228230503355e-5 ± 6.459567302786832e-6
The "SNR" is 12.025926608353366

Visualize

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