直交座標と楕円体高から緯度経度とジオイド高を計算する

6647 ワード
pyproj 緯度経度 地理院地図 地理院地図 テキストリンク

直交座標(X,Y)と楕円体高Hが与えられているとします。
まず、下の関数で直交座標(X,Y)から緯度経度(Lon,Lat)に変換します。
座標系は
http://tmizu23.hatenablog.com/entry/20091215/1260868350
https://www.gsi.go.jp/sokuchikijun/jpc.html
あたりを参考に。

convert.py
def Convert(X,Y):
    EPSG4612 = pyproj.Proj("+init=EPSG:4612") #GRS80楕円体（緯度経度）
    EPSG2450 = pyproj.Proj("+init=EPSG:2450") #日本平面直角座標系8系
    return pyproj.transform(EPSG2450, EPSG4612, X,Y)

さて、標高は楕円体高からジオイド高を引くことで得られるので、ジオイド高を求める必要があります。国土地理院のジオイド高計算用のAPIを使うこととします。
https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/calcframe.html
https://vldb.gsi.go.jp/sokuchi/surveycalc/api_help.html

geoid.py
def Getgeoid(latitude,longtitude):
    base_url='https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/cgi/geoidcalc.pl'
    outputType='json'
    url=base_url + '?outputType=' + outputType + '&latitude=' + str(latitude) + '&longitude=' + str(longtitude)
    for j in range(10):
        try:
            req = urllib.request.Request(url)
            body = json.load(urllib.request.urlopen(req))
            geoidHeight = float(body['OutputData']['geoidHeight'])
            return geoidHeight
        except :
            print("再試行")
            time.sleep(2)
            pass
        else:
            break
    else:
        print("再試行全失敗")
    return geoidHeight

if __name__ == '__main__':
    N= #計算数
    geoidHeight=np.zeros(N)
    for i in range(リスト数):
        geoidHeight[i] = Getgeoid(latitude[i],longtitude[i])
        time.sleep(1)

time.sleepは地理院サーバーへの思いやり