データ
references.json
座標系の決定に使う基準星のデータです。
[
{
"name": "BD-16 1589",
"x": 3826,
"y": 3398,
"ra": "06h44m58.68333s",
"dec": "-16d47m50.6410s"
},
{
"name": "BD-16 1586",
"x": 5149,
"y": 2276,
"ra": "06h44m39.08151s",
"dec": "-16d43m43.7368s"
},
{
"name": "Gaia DR2 2947057475918406784",
"x": 2457,
"y": 888,
"ra": "06h45m19.97874s",
"dec": "-16d38m53.3484s"
}
]
位置計算の対象となる星のデータです。"ra", "dec" (天球座標)が書いてある星は "x", "y" (写真上の位置)の計算対象、"x", "y" が書いてある星は "ra", "dec" の計算対象です。
座標系が正しく決定されているか確認するために天球座標のわかっているいくつかの星のデータを追加してあります。
[
{
"name": "Sirius A (2007)",
"ra" : "06h45m08.91728s",
"dec" : "-16d42m58.0171s"
},
{
"name": "Sirius A (Stellarium 0.22.1 (2022/3/12))",
"ra" : "06h45m08.45s",
"dec" : "-16d43m39.2s"
},
{
"name": "Sirius A (observed (2022/3/12))",
"x": 3240,
"y": 2161
},
{
"name": "Sirius B (observed (2022/3/12))",
"x": 3193,
"y": 2135
},
{
"name": "TYC 5949-2705-1",
"ra" : "06h45m39.85092s",
"dec" : "-16d37m27.3042s"
},
{
"name": "TYC 5949-2765-1",
"ra": "06h45m47.85991s",
"dec": "-16d41m08.7441s"
},
{
"name": "CSS 254",
"ra": "06h44m28.87639s",
"dec": "-16d42m58.3959s"
},
{
"name": "TYC 5949-2747-1",
"ra": "06h44m29.28203s",
"dec": "-16d37m20.2212s"
},
{
"name": "TYC 5949-2643-1",
"ra": "06h44m32.59814s",
"dec": "-16d35m32.3589s"
}
]
astro-image-masurement.py
references.json と targets.json と元画像を与えると計算結果(results.json)と星のマーカーとグリッドを記入したSVG画像を出力するスクリプトです。
import json
import math
import numpy as np
from astropy import units as u
from astropy.coordinates import Angle
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.wcs.utils import fit_wcs_from_points
import svgwrite
from PIL import Image
import base64
from argparse import ArgumentParser
argparser = ArgumentParser(description='get RA/DEC of target position.')
argparser.add_argument('ref_json', metavar='reference_stars.json',
help="reference stars data file.")
argparser.add_argument('target_json', metavar='target_stars.json',
help="target stars data file.")
argparser.add_argument('in_image', help="input jpg/png file.")
argparser.add_argument('out_svg', help="output svg file.")
args = argparser.parse_args()
ref_stars = None
target_stars = None
with open(args.ref_json, 'r', encoding='utf-8') as f:
ref_stars = json.load(f)
with open(args.target_json, 'r', encoding='utf-8') as f:
target_stars = json.load(f)
def ra_dec_to_angle(str_or_deg):
if (type(str_or_deg) is str):
return Angle(str_or_deg)
else:
return Angle(str_or_deg, u.deg)
ref_ra = []
ref_dec = []
ref_x = []
ref_y = []
proj_point = 'center'
for star in ref_stars:
ra = ra_dec_to_angle(star['ra'])
dec = ra_dec_to_angle(star['dec'])
ref_ra.append(ra.deg)
ref_dec.append(dec.deg)
ref_x.append(star['x'])
ref_y.append(star['y'])
if ('proj_point' in star) and star['proj_point']:
proj_point = SkyCoord(ra=ra.deg, dec=dec.deg, frame="icrs", unit=u.deg)
target_xy = []
xy_index = []
target_radec = []
radec_index = []
for i, star in enumerate(target_stars):
if 'x' in star:
xy_index.append(i)
target_xy.append([star['x'], star['y']])
elif 'ra' in star:
radec_index.append(i)
target_radec.append((ra_dec_to_angle(star['ra']).deg,
ra_dec_to_angle(star['dec']).deg))
stars = SkyCoord(ra=ref_ra, dec=ref_dec, frame="icrs", unit=u.deg)
pixels_x = np.array(ref_x)
pixels_y = np.array(ref_y)
wcs = fit_wcs_from_points(xy=[pixels_x, pixels_y],
world_coords=stars,
proj_point=proj_point)
radec_result = wcs.wcs_pix2world(target_xy, 0)
for i, star_index in enumerate(xy_index):
star = target_stars[star_index]
star['ra'] = Angle(radec_result[i][0], u.deg).to_string(unit=u.hour)
star['dec'] = Angle(radec_result[i][1], u.deg).to_string(unit=u.degree)
xy_result = wcs.wcs_world2pix(target_radec, 0)
for i, star_index in enumerate(radec_index):
star = target_stars[star_index]
star['x'] = xy_result[i][0]
star['y'] = xy_result[i][1]
output = json.dumps(target_stars, indent=2, ensure_ascii=False)
print(output)
image = Image.open(args.in_image)
w = image.width
h = image.height
corner_points = wcs.wcs_pix2world([[0, 0], [w, 0], [0, h], [w, h]], 0)
ra_min = Angle(min(list(map(lambda p: p[0], corner_points))), u.deg)
ra_max = Angle(max(list(map(lambda p: p[0], corner_points))), u.deg)
dec_min = Angle(min(list(map(lambda p: p[1], corner_points))), u.deg)
dec_max = Angle(max(list(map(lambda p: p[1], corner_points))), u.deg)
ra1 = ra_min.hms
ra2 = ra_max.hms
dec1 = dec_min.dms
dec2 = dec_max.dms
ra_h_start = int(ra1.h)
ra_h_last = int(ra2.h)
ra_m_start = int(ra1.m)
ra_m_last = int(ra2.m)
ra_s_start = int(ra1.s)
ra_s_last = int(ra2.s)
ra_scale = []
ra_scale_deg = []
for ra_h in range(ra_h_start, ra_h_last + 1):
if ra_h == ra_h_last:
ra_m_stop = ra_m_last + 1
else:
ra_m_stop = 60
for ra_m in range(ra_m_start, ra_m_stop):
if (ra_h == ra_h_last) and (ra_m == ra_m_last):
ra_s_stop = ra_s_last + 1
else:
ra_s_stop = 60
for ra_s in range(ra_s_start, ra_s_stop):
ra = (ra_h, ra_m, ra_s)
ra_scale.append(ra)
ra_scale_deg.append(Angle(ra, unit="hourangle").deg)
ra_s_start = 0
ra_m_start = 0
dec_d_start = int(dec1.d)
dec_d_last = int(dec2.d)
dec_m_start = int(dec1.m)
dec_m_last = int(dec2.m)
dec_s_start = int(dec1.s)
dec_s_last = int(dec2.s)
dec_scale = []
dec_scale_deg = []
for dec_d in range(dec_d_start, dec_d_last + 1):
if dec_d == dec_d_last:
dec_m_stop = dec_m_last + 1
else:
dec_m_stop = 60 if (dec_d >= 0) else 1
for dec_m in range(dec_m_start, dec_m_stop):
if (dec_d == dec_d_last) and (dec_m == dec_m_last):
dec_s_stop = dec_s_last + 1
else:
dec_s_stop = 60 if (dec_d >= 0) else 1
for dec_s in list(range(dec_s_start, dec_s_stop)):
dec = (dec_d, dec_m, dec_s)
dec_scale.append(dec)
dec_scale_deg.append(Angle(dec, unit=u.deg).deg)
dec_s_start = 0 if (dec_d >= 0) else -59
dec_m_start = 0 if (dec_d >= 0) else -59
drw = svgwrite.Drawing(args.out_svg, size=(w, h))
with open(args.in_image, 'rb') as f:
mime = ''
if args.in_image.upper().endswith('.JPG'):
mime = 'image/jpeg'
elif args.in_image.upper().endswith('.PNG'):
mime = 'image/png'
image_data = 'data:' + mime + ';base64,'
image_data += base64.standard_b64encode(f.read()).decode()
drw.add(drw.image(image_data))
ra_lines_h = []
ra_lines_m = []
ra_lines_s = []
for i,ra in enumerate(ra_scale):
coords = []
for dec_deg in dec_scale_deg:
coords.append((ra_scale_deg[i], dec_deg))
line = wcs.world_to_pixel(SkyCoord(coords, frame="icrs", unit=u.deg))
ra_h, ra_m, ra_s = ra
if (ra_m == 0) and (ra_s == 0):
ra_lines_h.append(line)
elif ra_s == 0:
ra_lines_m.append(line)
else:
ra_lines_s.append(line)
dec_lines_d = []
dec_lines_m = []
dec_lines_s = []
for i,dec in enumerate(dec_scale):
coords = []
for ra_deg in ra_scale_deg:
coords.append((ra_deg, dec_scale_deg[i]))
line = wcs.world_to_pixel(SkyCoord(coords, frame="icrs", unit=u.deg))
dec_d, dec_m, dec_s = dec
if (dec_m == 0) and (dec_s == 0):
dec_lines_d.append(line)
elif dec_s == 0:
dec_lines_m.append(line)
else:
dec_lines_s.append(line)
def draw_lines(drw, lines, style):
g = drw.g()
for line in lines:
points = []
for i, x in enumerate(line[0]):
points.append((x, line[1][i]))
svg_line = drw.polyline(points)
svg_line.update({ 'style' : style })
g.add(svg_line)
return g
grid_g = drw.g()
ra_g = drw.g()
ra_g.add(draw_lines(drw, ra_lines_h, "stroke: #fcc; stroke-width: 2;"))
ra_g.add(draw_lines(drw, ra_lines_m, "stroke: #ffc; stroke-width: 2;"))
ra_g.add(draw_lines(drw, ra_lines_s, "stroke: #ccc; stroke-width: 1;"))
grid_g.add(ra_g)
dec_g = drw.g()
dec_g.add(draw_lines(drw, dec_lines_d, "stroke: #fcc; stroke-width: 2;"))
dec_g.add(draw_lines(drw, dec_lines_m, "stroke: #ffc; stroke-width: 2;"))
dec_g.add(draw_lines(drw, dec_lines_s, "stroke: #ccc; stroke-width: 1;"))
grid_g.add(dec_g)
drw.add(grid_g)
def to_xy(pixel):
x = (pixel[0]).item()
y = (pixel[1]).item()
return (x, y)
def draw_markers(stars, color, font_size):
for star in stars:
x = star['x']
y = star['y']
l1 = drw.line(start=(x - 4, y), end=(x + 4, y),
stroke=color, stroke_width='1')
l2 = drw.line(start=(x, y - 4), end=(x, y + 4),
stroke=color, stroke_width='1')
drw.add(l1)
drw.add(l2)
t1 = drw.text(star['name'], x=[x+24], y=[y-(font_size * 2)],
fill=color, font_size=font_size)
t2 = drw.text(star['ra'], x=[x+24], y=[y-font_size],
fill=color, font_size=font_size)
t3 = drw.text(star['dec'], x=[x+24], y=[y],
fill=color, font_size=font_size)
g = drw.g()
g.add(t1)
g.add(t2)
g.add(t3)
drw.add(g)
draw_markers(ref_stars, "#8f8", 48)
draw_markers(target_stars, "#ff8", 24)
drw.save(pretty=True)
出力結果です。
[
{
"name": "Sirius A (2007)",
"ra": "06h45m08.91728s",
"dec": "-16d42m58.0171s",
"x": 3171.2872323943716,
"y": 2032.3461331144615
},
{
"name": "Sirius A (Stellarium 0.22.1 (2022/3/12))",
"ra": "06h45m08.45s",
"dec": "-16d43m39.2s",
"x": 3198.84000411331,
"y": 2223.511470771049
},
{
"name": "Sirius A (observed (2022/3/12))",
"x": 3240,
"y": 2161,
"ra": "6h45m07.84766535s",
"dec": "-16d43m25.56154398s"
},
{
"name": "Sirius B (observed (2022/3/12))",
"x": 3193,
"y": 2135,
"ra": "6h45m08.56224723s",
"dec": "-16d43m20.10826886s"
},
{
"name": "TYC 5949-2705-1",
"ra": "06h45m39.85092s",
"dec": "-16d37m27.3042s",
"x": 1143.6019365086354,
"y": 469.298876585322
},
{
"name": "TYC 5949-2765-1",
"ra": "06h45m47.85991s",
"dec": "-16d41m08.7441s",
"x": 593.4770269367737,
"y": 1486.5362595217844
},
{
"name": "CSS 254",
"ra": "06h44m28.87639s",
"dec": "-16d42m58.3959s",
"x": 5830.79180786849,
"y": 2077.298236424214
},
{
"name": "TYC 5949-2747-1",
"ra": "06h44m29.28203s",
"dec": "-16d37m20.2212s",
"x": 5833.743239682289,
"y": 511.1390967611835
},
{
"name": "TYC 5949-2643-1",
"ra": "06h44m32.59814s",
"dec": "-16d35m32.3589s",
"x": 5622.878222448987,
"y": 8.085561764319664
}
]
出力画像の方はこんな感じです。ただし、マーカーのテキストが重なる部分や、緑と黄色以外のマーカーの色は手動で(InkScapeで)修正してあります。
シリウスA+B (2022/3/12 19:52) 位置測定結果(緑の大きな文字のマーカーが基準点、その他の小さな文字のマーカーが
ピクセル位置からの計算結果)
入力画像は4秒露出画像の上に1/15秒露出画像を重ねたものです。1/15秒露出画像に微かに写った周囲の星を使って手動で位置合わせしています。
separation.py
計算結果(results.json)と、二つの星の "name" を指定して離角を計算するスクリプトです。
import sys
import json
from astropy import units as u
from astropy.coordinates import Angle
from astropy.coordinates import SkyCoord
from argparse import ArgumentParser
argparser = ArgumentParser(description='get RA/DEC of target position.')
argparser.add_argument('stars_json', metavar='stars.json',
help="stars data file.")
argparser.add_argument('star_name_a', help="name of star A.")
argparser.add_argument('star_name_b', help="name of star B.")
args = argparser.parse_args()
stars = None
with open(args.stars_json, 'r', encoding='utf-8') as f:
stars = json.load(f)
star_a = None
star_b = None
for star in stars:
if star['name'] == args.star_name_a:
star_a = star
elif star['name'] == args.star_name_b:
star_b = star
coord_a = SkyCoord(Angle(star_a['ra']).deg, Angle(star_a['dec']).deg, unit=u.deg)
coord_b = SkyCoord(Angle(star_b['ra']).deg, Angle(star_b['dec']).deg, unit=u.deg)
for_ra = SkyCoord(Angle(star_b['ra']).deg, Angle(star_a['dec']).deg, unit=u.deg)
for_dec = SkyCoord(Angle(star_a['ra']).deg, Angle(star_b['dec']).deg, unit=u.deg)
print(star_a['name'] + ": " + star_a['ra'] + " / " + star_a['dec'])
print(star_b['name'] + ": " + star_b['ra'] + " / " + star_b['dec'])
print("separation: " + coord_b.separation(coord_a).to_string(unit=u.deg) +
" (" + coord_a.separation(for_ra).to_string(unit=u.deg) + " / " +
coord_a.separation(for_dec).to_string(unit=u.deg) + ")")
results.json、"Sirius A (observed (2022/3/12))"、"Sirius B (observed (2022/3/12))" を指定し以下の出力を得ました。
Sirius A (observed (2022/3/12)): 6h45m07.84766535s / -16d43m25.56154398s
Sirius B (observed (2022/3/12)): 6h45m08.56224723s / -16d43m20.10826886s
separation: 0d00m11.62396968s (0d00m10.26536042s / 0d00m05.45327512s)
