J'essaie d'écrire un programme qui prendra un ensemble de longitude & coordonnées de latitude de l'utilisateur, les convertir en x & y coordonnées pour une carte de projection Mollweide, puis rapportez la valeur du pixel à ces coordonnées (dans ce cas, une température de bruit).Conversion latitude et longitude en x & y Mollweide coordonnées de la carte en utilisant une itération de Newton-Raphson en Python
La carte/données que j'utilise est la Haslam 408 MHz All Sky Survey qui est fournie sous forme de carte de projection Mollweide. Ces données sont au format .fits et constituent une vaste étude du bruit dans la bande 408 MHz.
Selon le Mollweide projection Wikipedia page, il est possible d'utiliser une itération Newton-Raphson pour convertir les coordonnées de longitude/latitude en coordonnées x/y. J'ai basé le schéma d'itération dans mon programme en grande partie sur les méthodes de la page de Wikipedia et dans this GitHub post. Toutefois, mon programme ne semble pas indiquer les valeurs correctes pour la longitude et la latitude que j'entre. Je pense surtout que l'un des deux (ou les deux) facteurs contribuent à cette erreur:
- La façon dont je suis en œuvre le système d'itération est incorrect, et entraînant ainsi des valeurs incorrectes signalées.
- Je ne comprends pas bien ce que la valeur de rayon, R, représente dans le schéma d'itération. Je ne trouve pas de littérature sur la façon de déterminer la valeur R appropriée au-delà de laquelle "R est le rayon du globe à projeter". J'ai supposé que cela serait basé sur la taille de la carte en pixels; dans ce cas, l'image de la carte est 4096x2048 pixels, donc j'ai essayé d'utiliser 2048, 1024, et simplement 1 comme les valeurs R, en vain.
Ci-dessous j'ai fourni mon code pour examen:
from math import sin, cos, pi, sqrt, asin
from astropy.io import fits
hdulist = fits.open('data.fits')
hdulist.info()
data = hdulist[1].data
sqrt2 = sqrt(2)
def solveNR(lat, epsilon=1e-6): #this solves the Newton Raphson iteration
if abs(lat) == pi/2:
return lat # avoid division by zero
theta = lat
while True:
nexttheta = theta - (
(2 * theta + sin(2 * theta) - pi * sin(lat))/
(2 + 2 * cos(2 * theta))
)
if abs(theta - nexttheta) < epsilon:
break
theta = nexttheta
return nexttheta
def checktheta(theta, lat): #this function is also currently unused while debugging
return (2 * theta + sin(2 * theta), pi * sin(lat))
def mollweide(lat, lon, lon_0=0, R=1024):
lat = lat * pi/180
lon = lon * pi/180
lon_0 = lon_0 * pi/180 # convert to radians
theta = solveNR(lat)
return (R * 2 * sqrt2 * (lon - lon_0) * cos(theta)/pi,
R * sqrt2 * sin(theta))
def inv_mollweide(x, y, lon_0=0, R=1024, degrees=True): # inverse procedure (x, y to lat, long). Currently unused
theta = asin(y/(R * sqrt2))
if degrees:
factor = 180/pi
else:
factor = 1
return (
asin((2 * theta + sin(2 * theta))/pi) * factor,
(lon_0 + pi * x/(2 * R * sqrt(2) * cos(theta))) * factor
)
def retrieve_temp(lat, long): #retrieves the noise temp from the data file after calling the mollweide function
lat = int(round(lat))
long = int(round(long))
coords = mollweide(lat, long)
x, y= coords
x = int(round(x))
y= int(round(y))
x = x-1
y = y-1
if x < 0:
x = x*(-1)
if y < 0:
y = y*(-1)
print("The noise temperature is: ",data[y, x],"K")
def prompt(): #this is the terminal UI
cont = 1
while cont == 1:
lat_cont = 1
while lat_cont == 1:
lat = float(input('Please enter the latitude: '))
lat_val = 1
while lat_val == 1:
if lat > 180 or lat < -180:
lat = float(input('Invalid input. Make sure your latitude value is in range -180 to 180 degrees \n'
'Please enter the latitude: '))
else:
lat_val = 0
lat_cont = 0
long_cont = 1
while long_cont == 1:
long = float(input('Please enter the longitude: '))
long_val = 1
while long_val == 1:
if long > 90 or long < -90:
long = float(input('Invalid input. Make sure your latitude value is in range -90 to 90 degrees \n'
'Please enter the latitude: '))
else:
long_val = 0
long_cont = 0
retrieve_temp(lat, long)
valid = 1
while valid == 1:
ans = input('Would you like to continue? Y or N: ').lower()
ans_val = 1
while ans_val ==1:
if not (ans == 'y' or ans == 'n'):
ans = input('Invalid input. Please answer Y or N to continue or exit: ')
elif ans == 'y':
ans_val = 0
cont = 1
valid = 0
elif ans == 'n':
ans_val = 0
cont = 0
valid = 0
prompt()
hdulist.close()
Toutes mes excuses si je ne suivre les conventions de Python typiques dans le code ci-dessus; Je suis nouveau à Python.
Juste pour confirmer, vous faites cela en tant que projet personnel, et vous ne voulez pas utiliser des toolkits existants, en particulier quelque chose comme proj4? – barrycarter
@barrycarter Oui, je suis un étudiant Astro de premier cycle qui fait cela en tant que projet d'été personnel pour garder mes compétences en codage avec quelque chose d'Astro. –