折角生成した wcs なので

さて、先日せっかく作った WCS なので具体的な遊びの一例を紹介したいと思います。 先日の WCS は 2023-10-23 23:00:00+09:00 に Atomcam2 で撮影した動画像から生成しました。 それを元に別の日時の画像のアノテーションを作成しようと思います。ただし、カメラの設置場所や向きを動かしておらず一定の撮影方向としていることを前提とします。

まず別日時の動画像を入手します。今回は 2023-10-27 00:15:00+09:00 の Atomcam2 の動画像を用います。

$ wget http://atomcam.local/sdcard/record/20231027/00/15.mp4

mp4 から加算画像を作成します。

1分程度加算しているので分かりづらいですが、やや曇り気味であまり星が見えていません。この画像ですと恐らくはプレートソルブも失敗しやすいと思います(試してないです)。

さて先日はここからプレートソルブをしたり、画面上の点と天球の座標のマップを作りアノテーションを作成したわけですが、今回は既に生成した WCS があるのでそれを流用します。 WCS には視点となる画像上の点の情報CRPIXとその点に対応する天球座標CRVALがあります。観測場所(緯度経度:LatLon)やカメラの向き(方位角仰角:AltAz)を一定とすると観測日時からその時点の新たな天球座標は計算可能になります。それを計算し WCS の CRVAL の値を更新し、その WCS を元にアノテーションを作成します。

from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import EarthLocation, AltAz, SkyCoord
import astropy.units as u

from datetime import datetime

# WCS 読み込み
fits = fits.open("new.wcs", "readonly")
wcs = WCS(header=fits[0].header)

# 元の WCS に対する観測日時を指定
obstime = datetime.fromisoformat("2023-10-23 23:00:30+09:00")
# 元の WCS に対する観測地点を指定(緯度経度)
location = EarthLocation(lon=1.397576692e02 * u.deg, 
                         lat=3.56802117e01 * u.deg, 
                         height=4 * u.m)

# WCS の CRVAL にある RaDec からカメラの視点としての AltAz を算出
altaz = AltAz(location=location, obstime=obstime)
view = SkyCoord(*wcs.wcs.crval, unit="deg").transform_to(altaz)

# 新しい観測日時を設定
obstime = datetime.fromisoformat("2023-10-27 00:15:30+09:00")
# 新しい観測日時を指定してカメラの視点である AltAz から RaDec を算出
radec = SkyCoord(az=view.az.deg, alt=view.alt.deg, 
                 unit="deg", frame="altaz", 
                 obstime=obstime, location=location).transform_to("icrs")

# 算出された RaDec を新たな CRVAL として WCS に上書き登録
fits[0].header["CRVAL1"] = radec.ra.deg
fits[0].header["CRVAL2"] = radec.dec.deg

# WCS をファイルに出力
fits.writeto('202310270015.wcs', overwrite=True)

出力した新たな WCS を用いてアノテーションを生成します

$ pngtopnm 202310270015.png | plot-constellations -w 202310270015.wcs -o test2.png -C -N -B -i -

うまくできているようです😀

atomcam をリニューアルしたので設置メモ/wcs の生成

暫く使っていた Atomcam2 が強烈な雨風に負けて水没してしまいリニューアルしました。

Atomcam2 購入

  • WiFi 接続初期設定、タイムスタンプ OFF
  • Atom ロゴ OFF
  • ナイトビジョンオート、LED Off (一度 on にして off)
  • 常時録画
  • 動画検出 OFF
  • sd カードフォーマット、atomcam tools 導入

適当な星空が撮影されるまで放置

撮影した動画像から wcs 生成
  • 星空が撮影された適当な動画像を Atomcam2 から取得
$ wget http://atomcam.local/sdcard/record/20231024/22/00.mp4
  • ノイズを軽減するために動画像から加算平均画像を生成
import cv2
import numpy as np

capture = cv2.VideoCapture("00.mp4")
width = int(capture.get(cv2.CAP_PROP_FRAME_WIDTH))
height = int(capture.get(cv2.CAP_PROP_FRAME_HEIGHT))
fps = int(capture.get(cv2.CAP_PROP_FPS))

r = np.zeros([height, width, 3], dtype=np.float32)
while True:
  (s, frame) = capture.read()
  if not s:
    break
  r += frame

r /= r.max()
cv2.imwrite("average.png", (r * 255).astype(np.uint8))

  • Astrometry.net のsolve-field をかけて画像の撮影範囲の目処をつけます
$ solve-field --scale-units arcsecperpix --scale-high 190 --scale-low 170 average.png  --overwrite --downsample 1

Atomcam2 のレンズは大体 180 arcsec/pixel 前後のようです。なので --scale-units で arcsecperpix を指定し、 --scale-high,--scale-low を用いて検索範囲を 170〜190 に限定することで solve-field での解決を早める期待ができます。 このコマンドラインで必ず解決できるわけではありません。特に歪みの少ない中心付近に解決に十分な星が検出できなければ解決に失敗する場合が多いと思います。 solve-field で無事解決ができると wcs と、アノテーションされた画像が取得できます。 これで問題なければ、以降出力された wcs を利用していけば良いのですが、残念ながら中心から離れるほどレンズの歪みの影響を受けて大きくずれていくという状態になるかと思います。 なのでここからはなるべく歪みに対応しズレが小さくなることが期待できる wcs を生成します。

プレートソルブせず人力で位置合わせをして wcs を生成

プレートソルブを用いて自動で位置を検出すると歪みに対応しきれないので手動で位置を設定して wcs を生成するようにします。

  • 画像に写っている星の画面上の (x, y) 座標位置を抽出

やり方は色々あるかと思いますが、ここでは安直に二値化して findContours して momtent で座標を決定しています。

import cv2
import pandas as pd

data = cv2.imread("average.png")
b = cv2.threshold(cv2.cvtColor(data, cv2.COLOR_BGR2GRAY), 127, 255, cv2.THRESH_BINARY)[1]
contours = cv2.findContours(b, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)[0]
stars = []
for c in contours:
    m = cv2.moments(c)
    if m['m00'] > 0:
        x,y= m['m10']/m['m00'] , m['m01']/m['m00']
        stars.append([x, y])
        cv2.circle(data, (round(x), round(y)), 2, (255, 0, 0), thickness=5)

pd.DataFrame(stars, columns=["x", "y"])

SEP とか使ったほうが良いのかな？

index x y
0 1893.6666666666665 1078.6666666666665
1 1676.3333333333333 1078.6666666666665
2 48.33333333333333 1078.6666666666665
3 37.5 1078.5
4 1916.9629629629626 1076.7037037037035
5 1866.3333333333333 1075.6666666666665
6 19.866666666666667 1075.8666666666666
7 1915.6538461538462 1072.4615384615386
8 769.6527777777777 1067.9166666666667
9 1533.9375 1052.3125
10 784.7619047619047 1039.5238095238094

後略

  • 抽出された星の(x, y) 座標に、その星自体の赤緯/赤経(ra, dec)の座標を紐付け fits 形式で保存します。 ここでは map.fits として保存しています。
from astropy.table import Table
import pandas as pd

map = pd.read_csv("map.csv", header=None, names=['INDEX_RA', 'INDEX_DEC', 'FIELD_X', 'FIELD_Y']).astype("double")
Table.from_pandas(map).write("map.fits", format="fits", overwrite=True)
index INDEX_RA INDEX_DEC FIELD_X FIELD_Y
0 89.93015897 37.21276409 769.2764227642276 1068.3495934959349
1 73.5628985 2.44067149 1534.048611111111 1052.125
2 87.7600428 37.30568175 784.6428571428571 1039.571428571429
3 81.57290804 28.60787346 989.452380952381 1030.2380952380952
4 72.80152507 5.60510146 1480.952380952381 1024.571428571429
5 87.87246934 39.14847936 751.3859649122807 1021.8070175438596
6 87.29358126 39.18113322 755.5555555555555 1014.2222222222222
7 72.45890935 6.96124744 1458.7395833333333 1011.375
8 72.65300943 8.90025258 1421.3333333333333 1005.1333333333332
9 96.22459158 49.28789903 518.6969696969696 995.6969696969696
10 89.88237261 44.94743492 636.4642857142857 982.1666666666666

後略、実際には 85 程のマップを作成しています。

  • Astrometry.net の fit-wcs を用いて wcs を生成します。
$ fit-wcs -W 1920 -H 1080 -s 4 -c map.fits -o new.wcs

-s オプションは SIP の歪み係数の次数の設定。どの程度が良いかはいろいろ試して出力を比較し決定しました。今回は概ね4を指定すると期待する出力に最も近くなりました。

  • 生成された wcs から Astrometry.net の plot-constellations を用いてアノテーションを生成します
$ pngtopnm average.png | plot-constellations -w new.wcs -o test.png -C -N -B -i -

こういった wcs で何ができるか

ずれの小さいアノテーションを被せられるようになり、頑張るとこういう動画が作れるようになったりします。

よしなしごと
  • 要は fit-wcs を使えばプレートソルブせず別途作成した座標のマップから wcs を生成できるという話なのだけれど、何故か fit-wcs に関するドキュメントが見つけられない。元々は SIP の対応を実装から読んでいてたまたま見つけて使ってみたらそれなりに動きそうという感じだった。