GPS经纬度坐标与XY坐标
一、Python源代码
def GPStoXY(self, lat, lon, ref_lat, ref_lon):
# input GPS and Reference GPS in degrees
# output XY in meters (m) X:North Y:East
lat_rad = math.radians(lat)
lon_rad = math.radians(lon)
ref_lat_rad = math.radians(ref_lat)
ref_lon_rad = math.radians(ref_lon)
sin_lat = math.sin(lat_rad)
cos_lat = math.cos(lat_rad)
ref_sin_lat = math.sin(ref_lat_rad)
ref_cos_lat = math.cos(ref_lat_rad)
cos_d_lon = math.cos(lon_rad - ref_lon_rad)
arg = np.clip(ref_sin_lat * sin_lat + ref_cos_lat * cos_lat * cos_d_lon, -1.0, 1.0)
c = math.acos(arg)
k = 1.0
if abs(c) > 0:
k = (c / math.sin(c))
x = float(k * (ref_cos_lat * sin_lat - ref_sin_lat * cos_lat * cos_d_lon) * self.CONSTANTS_RADIUS_OF_EARTH)
y = float(k * cos_lat * math.sin(lon_rad - ref_lon_rad) * self.CONSTANTS_RADIUS_OF_EARTH)
return x, y
二、公式推导
上面是将给定的Python程序转换为数学公式的过程：
设输入的GPS坐标为 (lat, lon)，参考GPS坐标为 (ref_lat, ref_lon)，输出的XY坐标为 (x, y)。
首先将经纬度转换为弧度制：
l
a
t
r
a
d
=
π
180
×
l
a
t
lat_{rad} = \frac{\pi}{180} \times lat
latrad=180π×lat
l
o
n
r
a
d
=
π
180
×
l
o
n
lon_{rad} = \frac{\pi}{180} \times lon
lonrad=180π×lon
r
e
f
l
a
t
r
a
d
=
π
180
×
r
e
f
l
a
t
ref_lat_{rad} = \frac{\pi}{180} \times ref_lat
reflatrad=180π×reflat
r
e
f
l
o
n
r
a
d
=
π
180
×
r
e
f
l
o
n
ref_lon_{rad} = \frac{\pi}{180} \times ref_lon
reflonrad=180π×reflon
接下来计算一些中间变量：
s
i
n
l
a
t
=
sin
(
l
a
t
r
a
d
)
sin_lat = \sin(lat_{rad})
sinlat=sin(latrad)
c
o
s
l
a
t
=
cos
(
l
a
t
r
a
d
)
cos_lat = \cos(lat_{rad})
coslat=cos(latrad)
r
e
f
s
i
n
l
a
t
=
sin
(
r
e
f
l
a
t
r
a
d
)
ref_sin_lat = \sin(ref_lat_{rad})
refsinlat=sin(reflatrad)
r
e
f
c
o
s
l
a
t
=
cos
(
r
e
f
l
a
t
r
a
d
)
ref_cos_lat = \cos(ref_lat_{rad})
refcoslat=cos(reflatrad)
KaTeX parse error: Double subscript at position 7: cos_d_̲lon = \cos(lon_…
然后根据以下公式计算XY坐标：
cos
(
c
)
=
sin
(
r
e
f
l
a
t
r
a
d
)
×
sin
(
l
a
t
r
a
d
)
+
cos
(
r
e
f
l
a
t
r
a
d
)
×
cos
(
l
a
t
r
a
d
)
×
cos
(
l
o
n
r
a
d
−
r
e
f
l
o
n
r
a
d
)
\cos(c) = \sin(ref_lat_{rad}) \times \sin(lat_{rad}) + \cos(ref_lat_{rad}) \times \cos(lat_{rad}) \times \cos(lon_{rad} - ref_lon_{rad})
cos(c)=sin(reflatrad)×sin(latrad)+cos(reflatrad)×cos(latrad)×cos(lonrad−reflonrad)
k
=
{
c
sin
(
c
)
if
∣
c
∣
>
0
1
if
∣
c
∣
=
0
k = \begin{cases} \frac{c}{\sin(c)} & \text{if } |c| > 0 \ 1 & \text{if } |c| = 0 \end{cases}
k={sin(c)cif ∣c∣>0 1if ∣c∣=0
x
=
k
×
(
r
e
f
c
o
s
l
a
t
×
sin
(
l
a
t
r
a
d
)
−
r
e
f
s
i
n
l
a
t
×
cos
(
l
a
t
r
a
d
)
×
cos
(
l
o
n
r
a
d
−
r
e
f
l
o
n
r
a
d
)
)
×
R
x = k \times (ref_cos_lat \times \sin(lat_{rad}) - ref_sin_lat \times \cos(lat_{rad}) \times \cos(lon_{rad} - ref_lon_{rad})) \times R
x=k×(refcoslat×sin(latrad)−refsinlat×cos(latrad)×cos(lonrad−reflonrad))×R
y
=
k
×
cos
(
l
a
t
r
a
d
)
×
sin
(
l
o
n
r
a
d
−
r
e
f
l
o
n
r
a
d
)
×
R
y = k \times \cos(lat_{rad}) \times \sin(lon_{rad} - ref_lon_{rad}) \times R
y=k×cos(latrad)×sin(lonrad−reflonrad)×R
其中 R 为地球半径，可以通过 self.CONSTANTS_RADIUS_OF_EARTH 获得。
因此，GPStoXY()函数的数学公式为：
KaTeX parse error: Double subscript at position 329: …_{rad}) \ cos_d_̲lon = \cos(lon_…
其中，输入的经纬度单位为度，输出的XY坐标单位为米。
扫一扫
1110
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
