1
我试图编写一个程序，该程序将从用户处获取一组经度为&的纬度坐标，将它们转换为x & y Mollweide投影映射的坐标，然后报告这些坐标处的像素值（在本例中为噪声温度）。将经度和纬度转换为x和y使用Python中的Newton-Raphson迭代转换Mollweide地图坐标
我正在使用的地图/数据是作为Mollweide投影图提供的Haslam 408 MHz All Sky Survey。这些数据采用.fits格式，是对408 MHz频段噪声的全天候调查。
根据Mollweide projection Wikipedia page，可以使用Newton-Raphson迭代将经度/纬度转换为x/y地图坐标。我在我的程序中基于维基百科页面和this GitHub post中的方法基于迭代方案。
但是，我的程序似乎没有报告我正在输入的经度和纬度的正确值。我基本上怀疑有两个（或两个）因素中的一个导致了这个错误：
- 我实现迭代方案的方式不正确，从而导致报告不正确的值。
- 我没有正确理解半径值R在迭代方案中的含义。我找不到任何关于如何确定超出“R是要投影的地球半径”的适当R值的文献。我认为这将基于像素大小的地图;在这种情况下，地图图像是4096x2048像素，所以我尝试过使用2048,1024，并且只是1作为R值，但无济于事。
下面我提供我的代码进行审核：如果
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()
道歉我没能按照上面的代码典型的Python公约;我是Python的新手。
只是为了确认，你是在做这个个人项目，而不想使用任何现有的工具包，特别是像proj4的东西？ – barrycarter
@barrycarter是的，我是一名本科Astro学生，将其作为一个个人夏季项目，以保持我的编码技能与Astro相关。 –