summaryrefslogtreecommitdiff
path: root/get_image.py
blob: e21e1e57f6bbe11ae64e7232ab9be339f83bf172 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#!/usr/bin/env python3
# encoding: utf-8
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi

"""
Get ugriz images from SDSS DR12 for given ra, dec combination in sexagesimal
format

Usage:

    python3 get_image.py [ra in degrees] [dec in degrees] [bands]

Excample:

    python3 get_image.py 150.0 45.0 u
    python3 get_image.py 150.0 45.0 ugriz
"""

import sys

import astropy.units as u
from astropy.coordinates import SkyCoord
from astroquery.sdss import SDSS
from astropy.wcs import WCS
import fitsio
import matplotlib as mpl
import matplotlib.pyplot as plt


def main():

    # Init
    ra = float(sys.argv[1])
    dec = float(sys.argv[2])
    band = sys.argv[3]

    # Query the region
    coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg)
    result = SDSS.query_region(coord,
                               radius=1 * u.arcmin,
                               data_release=12)

    # Fetch the first image of the query
    result = result[0]
    image = SDSS.get_images(run=result["run"],
                            camcol=result["camcol"],
                            field=result["field"],
                            data_release=12,
                            band=band)
    fitsio.write("tmp.fits", image[0][0])

    # Read in data, header, and a dx x dy sized cutout image. Note that the WCS
    # solution will only work on the original image.
    data = image[0][0].data
    header = image[0][0].header
    wcs = WCS(header)
    x, y = map(int, wcs.world_to_pixel(coord))
    dx = dy = 100
    cutout = data[y - dy // 2 : y + dy // 2,
                  x - dx // 2 : x + dx // 2]

    # Plot the first image, and a cutout of it as an example
    plt.figure(1)
    plt.imshow(data, origin="lower", norm=mpl.colors.LogNorm())
    plt.plot(x, y, 'rx')
    plt.gca().add_patch(
        plt.Rectangle([x - dx / 2, y - dy / 2], dx, dy, color="red", fc="none")
    )

    plt.figure(2)
    plt.imshow(cutout, norm=mpl.colors.LogNorm(), origin="lower")
    plt.show()


if __name__ == "__main__":
    main()