將經度和緯度轉換爲x和y使用Python中的Newton-Raphson迭代轉換Mollweide地圖座標

Question

我試圖編寫一個程序，該程序將從用戶處獲取一組經度爲&的緯度座標，將它們轉換爲x & y Mollweide投影映射的座標，然後報告這些座標處的像素值（在本例中爲噪聲溫度）。

我正在使用的地圖/數據是作爲Mollweide投影圖提供的Haslam 408 MHz All Sky Survey。這些數據採用.fits格式，是對408 MHz頻段噪聲的全天候調查。

根據Mollweide projection Wikipedia page，可以使用Newton-Raphson迭代將經度/緯度轉換爲x/y地圖座標。我在我的程序中基於維基百科頁面和this GitHub post中的方法基於迭代方案。

但是，我的程序似乎沒有報告我正在輸入的經度和緯度的正確值。我基本上懷疑有兩個（或兩個）因素中的一個導致了這個錯誤：

1. 我實現迭代方案的方式不正確，從而導致報告不正確的值。
2. 我沒有正確理解半徑值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的新手。

Answer 1

你的代碼看起來很合理。我的建議是要弄清楚什麼是錯誤的：

（1）嘗試評估你的mollweide和inv_mollweide函數，你知道結果應該是哪些點。例如。點赤道或主要子午線或類似的東西。

（2）您的mollweide和inv_mollweide實際上是反轉嗎？即如果你從你的輸出中取出你的輸出並把它放到另一箇中，你應該再次得到原始輸入。

（3）當您在地圖上移動時結果如何變化？你在某些地區（如地圖中部附近）獲得了正確的結果，而不是其他地區？當你接近邊緣會發生什麼？它是逐漸變得更加不準確還是有一些門檻，超過這個門檻你會得到非常不正確的答案？

我認爲牛頓方法的一個特點是，只有當你足夠接近解決方案才能收斂，否則你可以得到任何東西。我不知道你有多接近這個問題。

這似乎是一個很大的問題。祝好運並玩得開心點。