在旋转/缩放地图上平移纬度/经度坐标答案

【问题标题】：Translating lat/lon coordinate on rotated/scaled maps在旋转/缩放地图上平移纬度/经度坐标
【发布时间】：2013-06-16 14:43:27
【问题描述】：

我需要在旧地图中标出纬度/经度坐标，如下所示： http://www.davidrumsey.com/luna/servlet/detail/RUMSEY~8~1~3017~90020003:Topographical-Map-Of-The-City-and-C

在谷歌地图上投影： http://rumsey.geogarage.com/maps/g2784000.html

我知道地图的边界框以及图像的像素尺寸，但我不知道如何在此地图中定位纬度/经度。

左上角是 (40.698291,-74.079051)，左下角是 (40.659855,-73.979638)，右下角是 (40.855232,-73.835582)，右上角是 (40.882919,-73.940263)。

是否有关于如何在历史地图中定位地段/经度基准的标准公式？

【问题讨论】：

有四个角的坐标吗？
我用一个更好的例子更新了这个问题。酷吗？

标签： math maps

【解决方案1】：

只要您可以假设您的纬度/经度坐标形成一个矩形网格（这应该是城市规模的合理近似值，尽管对于较大的区域或非常精确的数据会失败），您可以假设您的转换在坐标为射影变换。请参阅this post（或this one，如果您更愿意继续使用SO）了解如何计算相应的变换矩阵并使用它来转换坐标的详细信息。

基于this post，我对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 和y 从零测量到一，其中零分别表示左侧。上边缘，一个在右边。下边缘。

如果矩形网格的假设还不够，您需要了解用于创建地图的投影的详细信息，以便将地图位置与球坐标相关联。

【讨论】：