经纬度坐标转换xy坐标 python_在Python中使用NewtonRaphson迭代将经纬度转换为xy Mo

我试图编写一个程序，从用户那里获取一组经度和纬度坐标，将它们转换为Mollweide投影图的x&y坐标，然后报告这些坐标处的像素值(在本例中，是噪声温度)。在

我使用的地图/数据是Haslam 408 MHz All Sky Survey，它是作为Mollweide投影图提供的。这是一个.fits格式的数据，是一个408MHz波段的大范围全天噪声调查。在

根据Mollweide projection Wikipedia page，可以使用Newton-Raphson迭代将经度/纬度转换为x/y地图坐标。我在程序中的迭代方案主要基于Wikipedia页面和this GitHub post中的方法。在

但是，我的程序似乎没有报告我输入的经度和纬度的正确值。我主要怀疑是两个(或两个)因素中的一个导致了这个错误：我实现迭代方案的方式是不正确的，因此导致报告了不正确的值。在

我不能正确理解迭代方案中半径值R代表什么。除了“R是要投影的地球的半径”之外，我找不到任何关于如何确定正确的R值的文献。我假设这将基于地图的像素大小；在本例中，地图图像是4096x2048像素，所以我尝试使用20481024和简单的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新手。在

经纬度坐标转换xy坐标 python_在Python中使用NewtonRaphson迭代将经纬度转换为xy Mollweide地图坐标...