使用python将纬度和经度转换为x和y网格系统答案

【问题标题】：convert latitude and longitude to x and y grid system using python使用python将纬度和经度转换为x和y网格系统
【发布时间】：2014-08-28 07:32:54
【问题描述】：

我有带有纬度和经度值的文件，我想将 x 和 y 转换为 km 我想测量到每个点的距离。

例如，我将经纬度的第一个点（分别为51.58，-124.6）

to (0,0) 在我的 x 和 y 系统中，所以基本上我想找出其他点是什么以及它们从原点开始的位置，所以我想找到 51.56(lat) -123.64(long) 是什么in (x,y) in km 等文件的其余部分。

我想在 python 中完成这一切，有一些排序代码吗？

例如，我在网上找到了网站

http://www.whoi.edu/marine/ndsf/cgi-bin/NDSFutility.cgi?form=0&from=LatLon&to=XY

确实想要我想做的，我只是不知道他们是怎么做到的。

【问题讨论】：

您确实意识到地球是一个球体，因此您的问题没有通用的解决方案（换句话说，如果没有一些近似/扭曲，世界无法映射到 2D 表面）。您可以在movable-type.co.uk/scripts/latlong.html 找到您需要的公式
如果你想尝试自己解决问题，你应该开始编码，当你有一个具体的、狭窄的问题时再回来。如果您只是想要一个现成的答案，这可能不是要查看的网站。
@Floris 在技术上是一个扁球体，但你的观点是正确的 :)
链接的网站是用HTML和js制作的，你可以把它翻译成python吗？
@Floris 在这种情况下，我假设与我的数据相对应的这个特定位置是平坦的，只是为了将 lat、long 转换为 x 和 y

标签： python latitude-longitude coordinate-systems

【解决方案1】：

您可以使用Great Circle Distance formula 获取 GPS 点之间的距离。纬度和经度位于大地坐标系中，因此您不能仅转换为平面 2D 网格并使用欧几里得距离。您可以将足够接近的点转换为近似网格，方法是采用像 (X,Y) 这样的任意点，将其设置为原点（就像您所做的那样），然后使用大圆距bearing 在平面上绘制相对于彼此的点，但这是一个近似值。

【讨论】：

【解决方案2】：

UTM 预测以米为单位。所以你可以在这个链接上使用类似 utm lib 的东西：

https://pypi.python.org/pypi/utm

用谷歌搜索 python lat lon 到 UTM 将指向几个选项。

UTM 区域的宽度为 6 度经度，从本初子午线的 0 开始。每个 UTM 带的原点位于赤道（x 轴）上，y 轴位于最西经度。这使得网格向北和向东偏正。您可以计算与这些结果的距离。 UTM 区域中间的值最准确。

您还应该知道原始纬度值所基于的数据，并在转换中使用相同的数据。

【讨论】：

【解决方案3】：

以下内容让您非常接近（以公里为单位回答）。如果你需要比这更好，你必须在数学上更加努力 - 例如通过遵循给出的一些链接。

import math
dx = (lon1-lon2)*40000*math.cos((lat1+lat2)*math.pi/360)/360
dy = (lat1-lat2)*40000/360

变量名应该很明显。这给了你

dx = 66.299 km (your link gives 66.577)
dy = 2.222 km (link gives 2.225)

选择坐标（例如，lon1, lat1）作为原点后，应该很容易了解如何计算所有其他 XY 坐标。

注意 - 系数 40,000 是地球的周长，以千米为单位（跨两极测量）。这让你接近。如果您查看您提供的链接的来源（您必须四处挖掘才能找到javascript which is in a separate file），您会发现他们使用了一个更复杂的公式：

function METERS_DEGLON(x)
{  
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return((111415.13 * cos(d2r))- (94.55 * cos(3.0*d2r)) + (0.12 * cos(5.0*d2r)));
   }
}

function METERS_DEGLAT(x)
{
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return(111132.09 - (566.05 * cos(2.0*d2r))+ (1.20 * cos(4.0*d2r)) - (0.002 * cos(6.0*d2r)));
   }
}

在我看来，他们实际上是在考虑地球并不完全是一个球体这一事实......但即便如此，当你做出假设时，你可以将地球的一部分视为你要去的平面有一些错误。我确信他们的公式错误更小......

【讨论】：

为什么需要乘以 40000
@learner 可能是因为地球的周长是 40075.017 公里（赤道）和 40007.86 公里（经向）
@yoshi 确实 - 我实际上在我的回答中这么说。我似乎记得有一次，米的定义是在两极测量的地球周长的 1/40000000。
我必须将 lon2-lon1 更改为 lon1-lon2 才能获得正确的结果。美国人会忽略经度的负数吗？
@R2-D2 看起来像是我的错字...显然dx 和dy 的顺序应该相同！

【解决方案4】：

如果您要使用 3D 系统，这些功能可以：

def arc_to_deg(arc):
    """convert spherical arc length [m] to great circle distance [deg]"""
    return float(arc)/6371/1000 * 180/math.pi

def deg_to_arc(deg):
    """convert great circle distance [deg] to spherical arc length [m]"""
    return float(deg)*6371*1000 * math.pi/180

def latlon_to_xyz(lat,lon):
    """Convert angluar to cartesian coordiantes

    latitude is the 90deg - zenith angle in range [-90;90]
    lonitude is the azimuthal angle in range [-180;180] 
    """
    r = 6371 # https://en.wikipedia.org/wiki/Earth_radius
    theta = math.pi/2 - math.radians(lat) 
    phi = math.radians(lon)
    x = r * math.sin(theta) * math.cos(phi) # bronstein (3.381a)
    y = r * math.sin(theta) * math.sin(phi)
    z = r * math.cos(theta)
    return [x,y,z]

def xyz_to_latlon (x,y,z):
    """Convert cartesian to angular lat/lon coordiantes"""
    r = math.sqrt(x**2 + y**2 + z**2)
    theta = math.asin(z/r) # https://stackoverflow.com/a/1185413/4933053
    phi = math.atan2(y,x)
    lat = math.degrees(theta)
    lon = math.degrees(phi)
    return [lat,lon]

【讨论】：

【解决方案5】：

您可以使用 UTM：

pip install utm

这是一个例子：

>>> import utm
>>> utm.from_latlon(51.2, 7.5)
(395201.3103811303, 5673135.241182375, 32, 'U')

返回的格式为(EASTING, NORTHING, ZONE_NUMBER, ZONE_LETTER)。

注意事项

它也适用于 NumPy 数组：

>>> utm.from_latlon(np.array([51.2, 49.0]), np.array([7.5, 8.4]))
(array([395201.31038113, 456114.59586214]),
 array([5673135.24118237, 5427629.20426126]),
 32,
 'U')

反过来：

>>> utm.to_latlon(340000, 5710000, 32, 'U')
(51.51852098408468, 6.693872395145327)

【讨论】：