2014-01-28 9 views
5

Ich habe gesucht und suchte nach einer Lösung für dieses Problem und bin dabei nichts.Kartesisches Projektionsproblem in einem FITS-Bild durch PyFITS/AstroPy

Ich erzeuge rechteckige FITS Bilder durch Matplotlib und wende anschließend WCS-Koordinaten an sie mit AstroPy (oder PyFITS). Meine Bilder befinden sich im galaktischen Breiten- und Längengrad. Daher sollten die für meine Karten geeigneten Header-Schlüsselwörter GLON-CAR und GLAT-CAR (für die kartesische Projektion) sein. Ich habe other maps betrachtet, die die gleiche Kartenprojektion in SAO DS9 verwenden, und die Koordinaten arbeiten groß ... das Gitter ist tadellos orthogonal, wie es sein sollte. Die FITS-Standardprojektionen können gefunden werden here.

Aber wenn ich meine Karten erzeuge, sind die Koordinaten überhaupt nicht kartesisch. Hier sehen Sie eine Gegenüberstellung meiner Karte (links) und einer anderen Referenzkarte der ungefähr gleichen Region (rechts). Beide sind GLON-CAR und GLAT-CAR in der FITS-Header aufgeführt, aber meins ist verrückt, wenn in SAO DS9 betrachtet (beachten Sie, dass das Koordinatengitter ist etwas, das SAO DS9 generiert basierend auf den Daten in der FITS-Header oder zumindest irgendwo in der FITS-Datei gespeichert)):

(left) my map, and (right) reference map. HEADER keywords are the same, both Cartesian

Dies ist problematisch, weil die koordinaten zuweisen Algorithmus falsche Koordinaten, die jedem pixel zugeordnet werden, wenn der Vorsprung nicht stimmt.

Hat jemand dies festgestellt, oder wissen, was könnte das Problem sein?

Ich habe versucht, andere Projektionen (nur um zu sehen, wie sie in SAO DS9 durchführen) und sie kommen gut aus ... aber meine Cartesian und Mercator Projektionen kommen nicht mit dem orthogonalen Gitter wie sie sollten.

Ich kann nicht glauben, dass dies ein Fehler in AstroPy wäre, aber ich kann keine andere Ursache finden ... außer meine Argumente in der Kopfzeile sind falsch formatiert, aber ich sehe immer noch nicht, wie das verursachen könnte das Problem, das ich erlebe. Oder würdest du etwas anderes empfehlen? (Ich habe mir die maplotlib-Grundkarte angesehen, hatte aber einige Probleme damit, das auf meinem Computer zu verwenden).

Mein Header-Code ist unten:

from __future__ import division 
import numpy as np 
from astropy.io import fits as pyfits # or use 'import pyfits, same thing' 

#(lots of code in between: defining variables and simple calculations... 
#probably not relevant) 

header['BSCALE'] = (1.00000, 'REAL = TAPE*BSCALE + BZERO') 
header['BZERO'] = (0.0) 
header['BUNIT'] = ('mag ', 'UNIT OF INTENSITY') 
header['BLANK'] = (-100.00, 'BLANK VALUE') 
header['CRVAL1'] = (glon_center, 'REF VALUE POINT DEGR') #FIRST COORDINATE OF THE CENTER 
header['CRPIX1'] = (center_x+0.5, 'REF POINT PIXEL LOCATION') ## REFERENCE X PIXEL 
header['CTYPE1'] = ('GLON-CAR', 'COORD TYPE : VALUE IS DEGR') 
header['CDELT1'] = (-glon_length/x_length, 'COORD VALUE INCREMENT WITH COUNT DGR') ### degrees per pixel 
header['CROTA1'] = (0, 'CCW ROTATION in DGR')    
header['CRVAL2'] = (glat_center, 'REF VALUE POINT DEGR') #Y COORDINATE OF THE CENTER 
header['CRPIX2'] = (center_y+0.5, 'REF POINT PIXEL LOCATION') #Y REFERENCE PIXEL 
header['CTYPE2'] = ('GLAT-CAR', 'COORD TYPE: VALUE IS DEGR') # WAS CAR OR TAN 
header['CDELT2'] = (glat_length/y_length, 'COORD VALUE INCREMENT WITH COUNT DGR') #degrees per pixel 
header['CROTA2'] = (rotation, 'CCW ROTATION IN DEGR')        #NEGATIVE ROTATES CCW around origin (bottom left). 
header['DATAMIN'] = (data_min, 'Minimum data value in the file') 
header['DATAMAX'] = (data_max, 'Maximum data value in the file') 
header['TELESCOP'] = ("Produced from 2MASS") 

pyfits.update(filename, map_data, header) 

Vielen Dank für jede Hilfe, die Sie anbieten können.

Antwort

7

In der modernen Definition der -CAR Projektion (von Calabretta et al.), Produziert GLON-CAR/GLAT-CAR Projektion nur ein geradliniges Gitter, wenn CRVAL2 auf Null gesetzt wird. Wenn CRVAL2 nicht Null ist, dann ist das Gitter gekrümmt (dies sollte nichts mit Astropy zu tun haben). Sie können versuchen, dies zu beheben, indem Sie CRVAL2 und CRPIX2 so anpassen, dass CRVAL2 Null ist. Hilft das?

Nur um zu klären, was ich meine, versuchen Sie, nachdem Sie den Code oben, und vor dem Schreiben der Datei aus:

header['CRPIX2'] -= header['CRVAL2']/header['CDELT2'] 
header['CRVAL2'] = 0. 

Glück gehabt?

Wenn Sie sich die Kopfzeile der Referenzdatei ansehen, die Sie angeschaut haben, sehen Sie, dass CRVAL2 dort Null ist. Nur um klar zu sein, es gibt nichts falsch mit CRVAL2 ist nicht Null, aber das Raster ist dann nicht mehr geradlinig.

+1

Das ist es !!!Wow, vielen Dank ... hätte diese Lösung nie alleine entdeckt. Nur ein Hinweis für jedermann auf diese in der Zukunft stolpern ... da dies bedeutet, dass die Referenz als Koordinatenursprung Einstellung (0,0), das letzte Stück des Puzzles ist den Pixelwert (CRPIX1 und CRPIX2) zu bestimmen, die für der Ursprung. Abhängig von den Koordinaten der Karte können also eines oder beide Referenzpixel negative Werte haben. Ich musste den Code von Astrofrog leicht modifizieren, damit er für mich funktionierte (insbesondere musste ich mit dem Inversen multiplizieren - also 1/header ['CDELT1'] - um die richtigen Werte zu erhalten. – Teachey

+1

@Tachey - ah ja, ich habe es um den falschen Weg gebracht - reparieren! – astrofrog