lgrcia / twirl Goto Github PK

View Code? Open in Web Editor NEW

44.0 3.0 7.0 2.03 MB

Astrometric plate solving package for Python.

Home Page: https://twirl.readthedocs.io

License: MIT License

Python 97.63% Dockerfile 2.31% Shell 0.06%

astronomy detector observation plate-solving python stars telescope

twirl's Introduction

twirl

Python package for astrometric plate solving

When the coordinates of an image (RA, dec) and the size of its field of view is approximately known, twirl can be used to compute a World Coordinate System (WCS) using GAIA reference stars.

Installation

twirl can be installed using pip:

pip install twirl

or using poetry:

poetry add twirl

Example Usage

twirl is designed to be complementary to the astropy package. It is used to compute a WCS by matching an image detected stars with catalog stars. To query the catalog stars, the center coordinates of the image must be known, as well as the size of the field of view. If not available, the latter can be computed using the known pixel scale of the detector.

Hence, the process starts by extracting the image RA-DEC center equatorial coordinate and compute the instrument field of view

import numpy as np

from astropy.io import fits
from astropy import units as u
from astropy.coordinates import SkyCoord

# Open some FITS image
hdu = fits.open("...")[0]

# get the center of the image
ra, dec = hdu.header["RA"], hdu.header["DEC"]
center = SkyCoord(ra, dec, unit=["deg", "deg"])

# and the size of its field of view
pixel = 0.66 * u.arcsec  # known pixel scale
shape = hdu.data.shape
fov = np.max(shape) * pixel.to(u.deg)

We can then query the gaia stars in the field using this information

import twirl

sky_coords = twirl.gaia_radecs(center, 1.2 * fov)[0:12]

and match the queried stars to stars detected in the image

# detect stars in the image
pixel_coords = twirl.find_peaks(hdu.data)[0:12]

# compute the World Coordinate System
wcs = twirl.compute_wcs(pixel_coords, sky_coords)

leading to a World Coordinate System object.

A more complete example is provided in docs/ipynb/wcs.ipynb

Development

Project Requirements

Python 3.11.*
Poetry for Python package and environment management.

Installing Dependencies

The twirl project manages Python package dependencies using Poetry. You'll need to follow the instructions for installation there.

Then you can start a shell session with the new environment with:

$ poetry shell

N.B. For development with vscode you will need to run the following command:

$ poetry config virtualenvs.in-project true

This will installed the poetry .venv in the root of the project and allow vscode to setup the environment correctly for development.

To start development, install all of the dependencies as:

$ poetry install

N.B. Ensure that any dependency changes are committed to source control, so everyone has a consistenct package dependecy list.

Acknowledgements

This package has made use of the algorithm from

Lang, D. et al. (2010). Astrometry.net: Blind Astrometric Calibration of Arbitrary Astronomical Images. The Astronomical Journal, 139(5), pp.1782–1800. doi:10.1088/0004-6256/139/5/1782.

implemented in

Garcia, L. J. et al. (2022). prose: a Python framework for modular astronomical images processing. MNRAS, vol. 509, no. 4, pp. 4817–4828, 2022. doi:10.1093/mnras/stab3113.

See this documentation page for the BibTeX entries.

twirl's People

Contributors

Stargazers

Watchers

Forkers

juliotux pierfra-rocci erschmny observerly xiaozangyu ysyy88 ppp-one

twirl's Issues

Fail to find transform

from #105

The sets of points:

import numpy as np

pixels = np.array([[1874.23682135,  906.47045569,  744.72233708,  862.92482627,
        1554.26186158,  377.22820759,   31.92397294,   66.8977319 ,
        1705.18939397,  386.0952344 , 1435.64933862, 1353.8251716 ,
        1534.62432011,  425.8644209 ,  386.17091984],
       [ 946.4496088 ,  460.92037899,  647.68930592,  867.8975847 ,
          22.10989145,  699.83604197, 1521.24951425, 1924.84590848,
         764.02973769,  396.82036351,  525.75953215, 1201.6207975 ,
         493.05218508, 1869.49881979,  114.5057253 ]])

radecs = np.array([[274.73148949, 274.86490329, 274.84561483, 274.85110499,
        274.90849086, 274.90721201, 274.75127985, 274.78820829,
        274.77117239, 274.79333541, 274.75682744, 274.87694318,
        274.76899521, 274.80959496, 274.92856186],
       [-68.14639379, -68.15990274, -68.16806503, -68.15018564,
        -68.15770553, -68.17102645, -68.15448174, -68.12906399,
        -68.16645804, -68.13535781, -68.13438334, -68.14672887,
        -68.16049091, -68.17644741, -68.12362587]])

have a known affine transform

T = np.array(
    [[-1.18158524e-04,  8.90493019e-07,  2.74952189e+02],
     [ 3.05602827e-07,  4.40285151e-05, -6.81886584e+01],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]]
     )

but find_transform(radecs.T, pixels.T) fails

Make citation cff file

[FEATURE]: add types support with py.typed mypy support across workspace

Feature Request

Hey @lgrcia I am now using this library quite extensively and there will be a hard requirement for typing information to ensure type safety for an existing Python application.

Are you happy for me to start working on this?

Motivation

A strongly-typed version of twirl will allow us to begin work on upstream type support for prose.

Alternatives

There are no alternatives, the library will need to be sanely typed.

Additional context

Typing is a really useful way to build robust software, I would strongly advise this for avoiding obvious foot guns.

3 points asterisms

Try solving using 3 points asterism. Many part of the 4 points can be used after the matching is done

ci: Implement ruff linting with Black code formatting

Implement ruff automated linting on save with the Black code formatter.

License

Checklist

I've searched the project's issues.

Question

The README contains an image saying this repository is licensed under the MIT license, but there is no LICENSE file. I am not sure whether it is legally required, so to be on the save side, could you confirm that this is the case?

I'm curious for more detail about the project's motivations.

Checklist

I've searched the project's issues.

Question

I've just come across your project and I'm impressed by the polish - it really looks like you've put a ton of effort into it this year!

It made me very curious to know more about your underlying motivation. Of course I read right at the top that your goal is to take an image and its center RA/Dec, and compute a WCS for that image. But it seems to me like you've taken a perhaps more thorough than necessary approach to do that and I'm very curious why. Or maybe this is actually the path of least resistance?

First off, maybe my understanding of WCS is wrong - I'm not much of an astrophotographer (yet) so I'm not steeped in this stuff. I thought a WCS was basically a transformation matrix that would give you the image's RA, Dec, and orientation towards celestial north. Perhaps it includes distortion? I think maybe I read that somewhere?

You cite Lang et al, so presumably you made an active decision not to simply reuse Astrometry.net - what were those reasons?

It's quite clear that this is not a lost in space solver - it looks like the primary limitation there is that you run online queries of the Gaia database that require internet access?

Anyway congratulations on a very neat and well put together package!

Additional context

None

Add tests

Thread to discuss writing more testings.

Added some tests in 3fade46 based on finding the correct transform between two sets of points (in pixels). So not based on WCS.

There is a case study in the new documentation with a real image. This could be used for real life testings (WCS).

Plate solving failing at poles

Bug Report

When attempting to plate solve images taken near or at the poles, plate solving fails.

(This was run using my local db of gaia/tmass since it was a lot faster than querying Gaia directly. But the code below uses Gaia directly and produces the same result.)

Speaking to @lgrcia, we probably need to do some spherical coordinates transformation in twirl to fix the issue at the poles.

How To Reproduce

This script generates a basic synthetic image at any ra/dec using Gaia, taken through a SPECULOOS like field of view.
It then uses twirl to compute the wcs of the generated image.
Finally, it calculates the angular separation from the center of the image to where twirl thinks the center of the image is.

For comparison, I successfully ran the synthetic images of two fields through https://nova.astrometry.net/:

Trappist-1 : https://nova.astrometry.net/user_images/8803416#grid
90 Dec, 180 RA: https://nova.astrometry.net/user_images/8803414#grid

import twirl

from astropy.coordinates import SkyCoord
from astropy.wcs import utils

import numpy as np
from astropy.wcs import WCS

plate_scale = 0.35  # [arcsec/pixel] known plate scale
numx = 2048 # size of image in x direction
numy = 2048 # size of image in y direction

def center_accuracy(ra, dec):
    '''
    Test the accuracy of twirl's WCS solution by comparing the it to the real center of a synthesised image.
    '''

    # center of image
    center = SkyCoord(ra=ra, dec=dec, unit="deg")

    # camera fov
    # fovx = (1/np.abs(np.cos(center.dec.rad)))*numx * plate_scale / 3600 # ra fov [deg]
    fovy = numy * plate_scale / 3600 # dec fov [deg]
    
    # create wcs object of simulated image
    wcs = WCS(naxis=2) # create wcs object for simulated image
    wcs.wcs.cdelt = [-plate_scale / 3600, -plate_scale / 3600]
    wcs.wcs.cunit = ["deg", "deg"]
    wcs.wcs.crpix = [int(numx/2), int(numy/2)]
    wcs.wcs.crval = [center.ra.deg, center.dec.deg]
    wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]

    # call gaia
    gaias = twirl.gaia_radecs(center, 1.5*fovy) # get gaia stars

    # convert gaia stars to pixel coordinates
    gaias_pixel = np.array(SkyCoord(gaias, unit="deg").to_pixel(wcs)).T

    gaias = gaias[0:18] # limit to 18 stars

    # calculate wcs of synthetic image
    stars = gaias_pixel[0:12] # limit to 12 stars
    wcs = twirl.compute_wcs(stars, gaias)
    wcs_center = utils.pixel_to_skycoord(numx/2, numy/2, wcs)

    return center.separation(wcs_center), gaias_pixel
    
def generate_image(ra, dec, output_file='test.fits'):

    from astropy.io import fits
    from scipy.ndimage import gaussian_filter

    ang_sep, gaias_pixel = center_accuracy(ra, dec)

    gaias_pixel = gaias_pixel.T

    # make base image with dark current and sky background
    base = np.zeros((numy, numx)) 

    # add stars to base, with fluxes as brightness
    for i in range(len(gaias_pixel[0])):
        try:
            base[int(gaias_pixel[1][i]), int(gaias_pixel[0][i])] += 2**16 - 1
        except:
            pass

    # convolve with gaussian to simulate seeing
    seeing = 2 # [arcsec]
    image = gaussian_filter(base, sigma=int((seeing/plate_scale)/2.355))

    # make image 16 bit
    image = image.astype(np.uint16)

    # save as fits
    hdu = fits.PrimaryHDU(image)
    hdr = hdu.header
    hdr['RA'] = ra
    hdr['DEC'] = dec
    hdu.writeto(output_file, overwrite=True)

if __name__ == "__main__":

    import json
    import matplotlib.pyplot as plt

    # generate images to compare with astrometry.net
    generate_image(ra=346.62083, dec=-5.04111, output_file="trappist-1.fits")
    generate_image(ra=180, dec=90, output_file="90dec_pole.fits")

    # test accuracy along a range of declinations, including the poles
    decs = np.linspace(-90, 90, 20)
    data = []
    for dec in decs:
        print(dec)
        ang, _ = center_accuracy(ra=180, dec=dec)
        data.append({"ra": 180, "dec": dec, "ang": ang.arcsec})

    # save data to file
    with open("data.json", "w") as f:
        json.dump(data, f)

    # plot results
    ra = [item['ra'] for item in data]
    dec = [item['dec'] for item in data]
    ang = [item['ang'] for item in data]

    fig = plt.figure(figsize=(16, 10))

    plt.scatter(dec, ang)

    plt.yscale('log')

    plt.xlabel('Dec [deg]')
    plt.ylabel('Angular Separation (real to platesolved center) [\"]')
    plt.title('WCS Error')

    plt.savefig('wcs_error.png')
    plt.show()

Environment

OS: macOS
Python version: 3.11

module 'twirl' has no attribute 'find_peaks'

Bug Report

I am getting this error

module 'twirl' has no attribute 'find_peaks'

that is not present with an earlier version of the code.

How To Reproduce

pip install twirl
python
import twirl
twirl.[TAB]

Code sample

Environment

OS: macOS

Screenshots

Expected behavior

Additional context

ValueError: attempt to get argmax of an empty sequence

In twirl.match.find_transform when no matches are found: ValueError: attempt to get argmax of an empty sequence

lgrcia / twirl Goto Github PK

twirl's Introduction

twirl

Installation

Example Usage

Development

Project Requirements

Installing Dependencies

Acknowledgements

twirl's People

Contributors

Stargazers

Watchers

Forkers

twirl's Issues

Feature Request

Motivation

Alternatives

Additional context

Checklist

Question

Checklist

Question

Additional context

Bug Report

How To Reproduce

Environment

Bug Report

How To Reproduce

Code sample

Environment

Screenshots

Expected behavior

Additional context

Recommend Projects

Recommend Topics

Recommend Org