radiofotograf
Aktives Mitglied
In einem anderen Thread hatte @Papa Brummbär gebeten, eine erste Hilfestellung zum Thema FITS-Dateien einlesen mit Python zu bekommen.
Für Leute, die noch gar keine Berührung mit Python hatten, stellt sich vermutlich als erstes die Frage, wie man Python selbst zum Laufen bekommt. Viele Leute installieren sich gar keine IDE, sondern arbeiten einfach mit einem normalen Texteditor und starten dann die entsprechenden Skripte auf Kommandozeile. Python hat aber auch eine sehr mächtige interaktive Shell, in die man gelangt, wenn man Python startet. Noch besser ist sogar IPython, weil es Syntaxvervollständigung hat und vieles mehr. Alternativ kann man sich auch einen sog. Jupyter Server hochziehen und dann bequem im Webbrowser "Notebooks" erstellen - ein Gemisch aus Programmcode, Bildern und Markdown Text. Aber es gibt auch klassische IDEs, z.B. Spyder (einsteigerfreundlich) oder PyCharm (für Profis).
In den meisten Fällen muss man Python jedenfalls erstmal installieren. Dazu würde ich gerade Anfängern raten, die Anaconda Python Distribution zu installieren. Die gibt es auch für alle Plattformen. Unter Linux ist Python meist schon installiert. Man sollte aber nicht die Systeminstallation nutzen, bzw. wenn dann neue Pakete (Python-Sprech für Bibliotheken) nur in den User-Space installieren. Am besten arbeitet man ohnehin mit virtuellen Pythonumgebungen, sobald es etwas professioneller werden soll, aber das führt hier etwas zu weit.
Das Gute an Anaconda ist, dass schon alles Nötige mitgeliefert wird: NumPy, SciPy, matplotlib - diese 3 Pakete bezeichnen den sog. Scientific Stack - aber auch Jupyter, IPython, Spyder, AstroPy (braucht man für FITS) u.v.m.
Unter Windows sollte man dann die Anaconda Shell starten. (Unter Linux reicht vermutlich irgendeine Shell.) Mit "python" oder "ipython" lässt sie sich starten. Mit Spyder kenne ich mich zu wenig aus, aber das sollte alternativ auch funktionieren.
In Python lädt man dann zunächst die nötigen Module:
Dann laden wir eine FITS-Datei ein. FITS-Dateien bestehen aus einer oder mehreren Header-Data-Units (HDU). Die Bilddateien, die in der Astrofotografie genutzt werden, bestehen fast immer nur aus einer HDU (der sog. Primary HDU). Außerdem gibt es verschiedene Typen: Images und Tabellen (Tabellen können ASCII oder Binärformat haben - aber das heißt nicht, dass die FITS-Datei selbst als ASCII abgespeichert ist!). Ihr werdet es schon vermutet haben, dass es sich bei Astrobildern um den "Image"-Typ handelt...
Folgendermaßen kann man das mit dem astropy FITS Modul lesen:
Im Prinzip war es das schon. Alles weitere ist dann nicht mehr FITS-spezifisch. Aber um wenigstens kurz zu zeigen, wie es weiter gehen könnte, hier noch ein paar Tips. Das eingelesene Header Objekt ist ein normales Python Dictionary und kann sehr leicht verwendet werden:
Das Datenobjekt ist ein NumPy Array, was sehr mächtige und schnelle Funktionen erlaubt. Das kann ich hier nicht im Detail erläutern, aber man kann z.B. sehr einfach statische Kenngrößen damit berechnen:
Je nachdem, wie das Array abgespeichert ist, hat man 2D oder 3D Daten. 2D ist es, wenn man Mono-Daten hat, oder noch nicht De-Bayered wurde. 3D kann vorkommen, wenn die RGB Kanäle in einzelnen Ebenen abgespeichert sind. In letzterem Fall möchte man dann meist den Mittelwert (oder was auch immer) pro Kanal haben, dann gibt man das sog. "Axis" Keyword mit an (dabei muss man aber vorher sicherstellen, dass die Dimensionen auch in der erwarteten Reihenfolge vorkommen, meist hat die erste Achse dann Länge 3 für die R, G und B Bilder:
Man kann auch schöne Plots mit matplotlib erstellen. Hier nur ganz simpel das Graustufen-Bild dargestellt:
Falls noch nicht De-Bayered wurde, geht das auch sehr fix, z.B.
Hier wurde die sehr mächtige NumPy Array Slice Technik angewendet.
Richtig toll ist auch Astropy. Man kann nicht nur Koordinaten transformieren sondern auch die Bildprojektion (im FITS/WCS Standard, WCS steht für World Coordinate System) automatisiert auswerten und z.B. Pixel-Koordinaten in physikalische Koordinaten umwandeln:
Für Leute, die noch gar keine Berührung mit Python hatten, stellt sich vermutlich als erstes die Frage, wie man Python selbst zum Laufen bekommt. Viele Leute installieren sich gar keine IDE, sondern arbeiten einfach mit einem normalen Texteditor und starten dann die entsprechenden Skripte auf Kommandozeile. Python hat aber auch eine sehr mächtige interaktive Shell, in die man gelangt, wenn man Python startet. Noch besser ist sogar IPython, weil es Syntaxvervollständigung hat und vieles mehr. Alternativ kann man sich auch einen sog. Jupyter Server hochziehen und dann bequem im Webbrowser "Notebooks" erstellen - ein Gemisch aus Programmcode, Bildern und Markdown Text. Aber es gibt auch klassische IDEs, z.B. Spyder (einsteigerfreundlich) oder PyCharm (für Profis).
In den meisten Fällen muss man Python jedenfalls erstmal installieren. Dazu würde ich gerade Anfängern raten, die Anaconda Python Distribution zu installieren. Die gibt es auch für alle Plattformen. Unter Linux ist Python meist schon installiert. Man sollte aber nicht die Systeminstallation nutzen, bzw. wenn dann neue Pakete (Python-Sprech für Bibliotheken) nur in den User-Space installieren. Am besten arbeitet man ohnehin mit virtuellen Pythonumgebungen, sobald es etwas professioneller werden soll, aber das führt hier etwas zu weit.
Das Gute an Anaconda ist, dass schon alles Nötige mitgeliefert wird: NumPy, SciPy, matplotlib - diese 3 Pakete bezeichnen den sog. Scientific Stack - aber auch Jupyter, IPython, Spyder, AstroPy (braucht man für FITS) u.v.m.
Unter Windows sollte man dann die Anaconda Shell starten. (Unter Linux reicht vermutlich irgendeine Shell.) Mit "python" oder "ipython" lässt sie sich starten. Mit Spyder kenne ich mich zu wenig aus, aber das sollte alternativ auch funktionieren.
In Python lädt man dann zunächst die nötigen Module:
Python:
import astropy.io.fits as pf
import numpy as np
import matplotlib.pyplot as pltDann laden wir eine FITS-Datei ein. FITS-Dateien bestehen aus einer oder mehreren Header-Data-Units (HDU). Die Bilddateien, die in der Astrofotografie genutzt werden, bestehen fast immer nur aus einer HDU (der sog. Primary HDU). Außerdem gibt es verschiedene Typen: Images und Tabellen (Tabellen können ASCII oder Binärformat haben - aber das heißt nicht, dass die FITS-Datei selbst als ASCII abgespeichert ist!). Ihr werdet es schon vermutet haben, dass es sich bei Astrobildern um den "Image"-Typ handelt...
Folgendermaßen kann man das mit dem astropy FITS Modul lesen:
Python:
hdus = pf.open('<PFAD zur FITS-Datei>')
phdu = hdus[0] # read Primary HDU
# Display header information:
print(phdu.header)
# Display data block:
print(phdu.data)Im Prinzip war es das schon. Alles weitere ist dann nicht mehr FITS-spezifisch. Aber um wenigstens kurz zu zeigen, wie es weiter gehen könnte, hier noch ein paar Tips. Das eingelesene Header Objekt ist ein normales Python Dictionary und kann sehr leicht verwendet werden:
Python:
# Display and work with variables from the header:
date_obs = phdu.header['DATE-OBS']
print(date_obs)
# Can convert this to astropy Time (useful for coordinate convertions):
from astropy.time import Time
date_obs_ap = Time(date_obs)
print(date_obs_ap, date_obs_ap.mjd) # conversion to MJD is easy
time_since_observation = Time.now() - date_obs_ap
print(time_since_observation.sec, " s")Das Datenobjekt ist ein NumPy Array, was sehr mächtige und schnelle Funktionen erlaubt. Das kann ich hier nicht im Detail erläutern, aber man kann z.B. sehr einfach statische Kenngrößen damit berechnen:
Python:
# Display mean, median and RMS of data:
print(np.mean(phdu.data))
print(np.median(phdu.data))
print(np.std(phdu.data))Je nachdem, wie das Array abgespeichert ist, hat man 2D oder 3D Daten. 2D ist es, wenn man Mono-Daten hat, oder noch nicht De-Bayered wurde. 3D kann vorkommen, wenn die RGB Kanäle in einzelnen Ebenen abgespeichert sind. In letzterem Fall möchte man dann meist den Mittelwert (oder was auch immer) pro Kanal haben, dann gibt man das sog. "Axis" Keyword mit an (dabei muss man aber vorher sicherstellen, dass die Dimensionen auch in der erwarteten Reihenfolge vorkommen, meist hat die erste Achse dann Länge 3 für die R, G und B Bilder:
Python:
# Display mean, median and RMS of 3D array, assuming first axis has length 3
print(np.mean(phdu.data), axis=0)
print(np.median(phdu.data), axis=0)
print(np.std(phdu.data), axis=0)Man kann auch schöne Plots mit matplotlib erstellen. Hier nur ganz simpel das Graustufen-Bild dargestellt:
Python:
data = phdu.data
plt.close()
plt.imshow(data, origin='lower', interpolation='nearest')
plt.show()Falls noch nicht De-Bayered wurde, geht das auch sehr fix, z.B.
Python:
channel1 = data[0::2, 0::2]
channel2 = data[1::2, 0::2]
channel3 = data[0::2, 1::2]
channel4 = data[1::2, 1::2]
# TODO: figure out which channel is which color...
plt.close()
plt.imshow(channel1, origin='lower', interpolation='nearest')
plt.show()Richtig toll ist auch Astropy. Man kann nicht nur Koordinaten transformieren sondern auch die Bildprojektion (im FITS/WCS Standard, WCS steht für World Coordinate System) automatisiert auswerten und z.B. Pixel-Koordinaten in physikalische Koordinaten umwandeln:
Python:
from astropy.wcs import WCS
my_wcs = WCS(phdu.header)
# FITS indexes are 1-based (as in Fortran); As NumPy Arrays are
# zero-based, we need to provide origin=0:
xpix, ypix = 20, 30
lon, lat = my_wcs.all_pix2world(xpix, ypix, 0)
print(lon, lat)
# reverse direction
lon, lat = 130.3, -1.
xpix, ypix = my_wcs.all_world2pix(lon, lat, 0)
print(xpix, ypix)
Zuletzt bearbeitet:
