AGO Teleskop

Pri obdelavah FITS slik (datotek) lahko uporabimo zelo močno programsko orodje python in modul astropy. Kako to naredimo je prikazano v primeru izračuna centroida zvezde. Program je tudi v mapi za izmenjavo.

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# to je glava python programa in nastavljen utf-8, da lahko tipkam čžš
# tudi naredim izvršljiv program z ukazom v bash konzoli:
# chmod a+x mojfits.py 
#
#
# Program za branje FITS datotek in izračun centroida
#
# Bojan Dintinjana 
# nov. 2018


import math
from astropy.io import fits

# odprem mojo fits sliko:
hdulist = fits.open('Caecilia_V_300.fts')

# izpis informacije o sliki
hdulist.info()


# branje headerja in izpis časa posnetka
# slike Vega imajo samo en HDU (Header Data Unit) 
# zato indeks [0]

hdr = hdulist[0].header
print hdr['DATE-OBS']


#zdaj branje podatkov v 2D array
data = hdulist[0].data


# približne začetne koordinate moje zvezde za meritev, ki jih dobim recimo v ds9, 
# kjer gledam koordinati physical X, Y

x,y=903,1021

# izberenm si kvadratek 15x15 okoli zvezde, lahko bi dal tudi manj
# tudi upoštevam, da ima prvi piksel v ds9 ali v IRAFu koordinate (1,1) to je levo spodaj
# in upoštevam, da ima python prvi piksel data[0,0] in običajno levo zgoraj
# izračunam meje kvadratka; 
x1 = x - 7
x2 = x + 8
y1 = y - 7 
y2 = y + 8 

# izpis za kontrolo tako, zapis je [y,x] koordinati sta zamenjani, 
# x hitrejša os, ki je v fits glavi NAXIS1 je v pythonu na drugem mestu
#
#data[1020-5:1020+6,903-5:903+6]

# izračunam centroid zvezde
# najprej nastavim vsote na 0.
M = 0. 
Sx= 0.
Sy= 0.

# zanka teče po y in x in računam
for j in range(y1, y2):
	for i in range(x1, x2):
		# lep izpis matrike za kontrolo
		print ("%5d "% data[j,i]),
		M  += data[j,i]
		Sx += data[j,i] * i
		Sy += data[j,i] * j
	print

# končni izračun
Rx = 1./M * Sx
Ry = 1./M * Sy

# in izpis		
print "M =%7.0f"%M
print "Sx=%7.0f Sy=%7.0f"%(Sx,Sy)
print "Rx=%7.3f Ry=%7.3f"%(Rx,Ry)

$ ./fits.py 
Filename: Caecilia_V_300.fts
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU     107   (2049, 2049)   int16 (rescales to uint16)   
2015-11-16T23:54:50.12
  502    539    586    579    587    588    670    610    661    652    587    511    580    515    518 
  549    551    596    576    650    750    780    840    753    705    668    597    513    544    524 
  546    609    637    698    961   1201   1262   1244   1151    927    757    627    573    528    537 
  581    679    786   1003   1428   1889   2226   2222   1904   1457   1009    742    611    646    565 
  588    648    875   1309   1978   3013   3907   4066   3296   2234   1337    916    683    611    526 
  640    806    997   1428   2456   3846   5409   6120   5075   3369   1865   1109    758    632    562 
  684    843   1032   1534   2553   3978   5870   6938   6342   4312   2447   1290    836    619    578 
  695    824   1074   1591   2237   3283   4673   6008   5959   4325   2412   1375    830    635    575 
  736    809   1116   1422   1873   2465   3199   3897   3949   3199   2008   1164    757    632    583 
  696    820   1028   1364   1578   1866   2058   2274   2371   2045   1302    933    691    654    533 
  659    755    912   1117   1321   1418   1450   1380   1397   1219    975    764    635    615    568 
  631    668    798    861    973   1085    966   1018    988    831    725    677    578    569    512 
  592    584    689    720    760    797    753    795    706    658    627    611    513    535    532 
  542    540    606    606    617    678    649    626    609    601    561    508    489    512    523 
  550    496    539    525    574    586    566    588    590    531    532    555    503    532    513 
M = 287504
Sx=259567400 Sy=293422902
Rx=902.831 Ry=1020.587

 903.85267 1021.5628

Kar se lepo ujema z izračunanim, upoštevati je potrebno še +1 piksel zaradi izhodišča, ki je v pythonu [0,0].