1

次のような古い地図で緯度/経度の座標を示す必要があります: http://www.davidrumsey.com/luna/servlet/detail/RUMSEY~8~1~3017~90020003:Topographical-Map-Of-The-シティアンドシー

Google マップに投影: http://rumsey.geogarage.com/maps/g2784000.html

マップの境界ボックスと画像のピクセル単位の寸法はわかっていますが、このマップ内で緯度/経度を特定する方法がわかりませんでした。

左上 (40.698291,-74.079051)、左下 (40.659855,-73.979638)、右下 (40.855232,-73.835582)、右上 (40.882919,-73.940263) です。

歴史的な地図でロット/経度データムを見つける方法に関する標準的な公式はありますか?

4

1 に答える 1

2

緯度/経度座標が長方形のグリッドを形成すると想定できる限り (これは、都市規模での妥当な概算ですが、より広いエリアや非常に正確なデータでは失敗します)、座標間の変換は次のようになると想定できます。射影変換。対応する変換行列を計算し、それを使用して座標を変換する方法の詳細については、この投稿 (SO にとどまりたい場合はこの投稿) を参照してください。

この投稿に基づいて、私はsageで少し計算を行いました:

def pm1(a, b, c, d):
    M = Matrix([a , b, c]).transpose()
    f = M.adjoint() * d
    return M*diagonal_matrix(f)
def pm(a1, a2, b1, b2, c1, c2, d1, d2):
    return pm1(a2, b2, c2, d2)*(pm1(a1, b1, c1, d1).adjoint())
P1 = vector(QQ, [0, 0, 1])
P2 = vector(QQ, [0, 1, 1])
P3 = vector(QQ, [1, 1, 1])
P4 = vector(QQ, [1, 0, 1])
Q1 = vector(QQ, [40698291, -74079051, 1000000]) # Upper left
Q2 = vector(QQ, [40659855, -73979638, 1000000]) # Lower left
Q3 = vector(QQ, [40855232, -73835582, 1000000]) # Lower right
Q4 = vector(QQ, [40882919, -73940263, 1000000]) # Upper right
M = pm(Q1, P1, Q2, P2, Q3, P3, Q4, P4)
M.change_ring(RDF)/1e40

結果の式は次のようになります。

z =   559.910562534*lat - 539.510073656*lon - 64123.7576703
x = (-5629.59680416*lat - 2176.56828347*lon + 67876.8560722)/z
y = ( 8466.61769096*lat - 11263.0392472*lon - 1178932.12938)/z

結果の座標はx、y0 から 1 まで測定されます。ここで、0 は左の resp を示します。上端、および右の1つ。下端。

長方形のグリッドの仮定が十分でない場合は、マップの作成に使用された投影法の詳細を知る必要があります。これにより、マップの位置を球面座標に関連付けることができます。

于 2013-06-17T13:05:16.697 に答える