如何在 Pandas 中将 NAD83 坐标转换为常规 GPS 坐标（纬度/经度）？

Question

我在Python / Pandas中找不到任何库或函数可以将坐标从NAD83标准转换为GPS纬度/经度。

以下是一个重现此情况的示例。我正在寻找一种将这些NAD83坐标转换为常规GPS纬度/经度的功能。下面所有点的期望位置都在圣路易斯附近（即大约在38.6270° N，90.1994° W）。

xcoords = [894557.5, 880625.2, 896551.8, 896551.8, 896497.6, 903061.1]
ycoords = [1025952, 996012.7, 1025333, 1025333, 997157.3, 1033547]
df = pd.DataFrame({'xcoords':xcoords, 'ycoords':ycoords}) 

这意味着

xcoords     ycoords 
894557.5    1025952.0
880625.2    996012.7
896551.8    1025333.0
896551.8    1025333.0
896497.6    997157.3
903061.1    1033547.0

数据源 = http://www.slmpd.org/Crimereports.shtml

元数据 = http://www.slmpd.org/Crime/CrimeDataFrequentlyAskedQuestions.pdf （请查看 XCoord 和 YCoord - 这些坐标是关于什么的？）

- Rotail

Answer 1

有点晚了，但是......以下是我如何将SLMPD的x,y坐标转换为经纬度的方法。首先需要确定我们所在的EPSG区域，这是“epsg 3602 NAD83(NSRS2007) / Missouri East”。可以从https://spatialreference.org/ref/epsg/获取。通过一般研究，我得出将其转换为espg：4326的结论，“3D系统的水平分量。GPS卫星导航系统和北约军事测量使用。”

from pyproj import Transformer 

#transformer to convert from (epsg 3602 NAD83(NSRS2007) / Missouri East) to 4326 (WGS84 - World Geodetic System 1984, used in GPS)

transformer = Transformer.from_crs( "epsg:3602","epsg:4326",always_xy=False)

#empty lists to hold converted x y coordinates
lat = []
lon = []

#loop through columns 'XCoord' and 'YCoord' and declare variables x and y to use in transformer. 
#note I had to change the column values from feet to meters for the transformer to correctly work
for index, row in df.iterrows():
    x = (row['XCoord'] / 3.28)
    y = (row['YCoord'] / 3.28)
    
    #call transformer
    x1, y1 = transformer.transform(x,y)
      
    #append to lists    
    lat.append(x1)
    lon.append(y1)
    
# add new columns to df     
df['Lat'] = lat
df['Lon'] = lon

- SpacePatroller

Answer 2

首先，我们需要确定圣路易斯市的SPCS或FIPS区域。从NOAA手册NPS NGS 5中，我们得到2401号区域。然后，我们从特殊参考下载相应的proj字符串（或者您可以从NOAA手册附录A第68页获取必要数据）。有了这些数据，我们可以使用pyproj将SPCS数据转换为GPS经纬度值：

import pandas as pd
import pyproj

xcoords = [894557.5, 880625.2, 896551.8, 896551.8, 896497.6, 903061.1]
ycoords = [1025952, 996012.7, 1025333, 1025333, 997157.3, 1033547]
df = pd.DataFrame({'xcoords':xcoords, 'ycoords':ycoords})

fips2401 = pyproj.Proj("+proj=tmerc +lat_0=35.83333333333334 +lon_0=-90.5 +k=0.9999333333333333 +x_0=250000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")
wgs84 = pyproj.Proj("+init=EPSG:4326")

df[['lon', 'lat']] = pd.DataFrame(pyproj.transform(fips2401, wgs84, df.xcoords.to_numpy(), df.ycoords.to_numpy())).T

结果：

    xcoords    ycoords        lon        lat
0  894557.5  1025952.0 -90.239655  38.650896
1  880625.2   996012.7 -90.288682  38.568784
2  896551.8  1025333.0 -90.232678  38.649181
3  896551.8  1025333.0 -90.232678  38.649181
4  896497.6   997157.3 -90.233154  38.571814
5  903061.1  1033547.0 -90.209794  38.671681

- Stef

---

我认为lon需要在+90左右而不是-90。这与现实相差一个负号！ - Rotail

哈！这真的让我感到困惑，因为现在运行它时我得到了不同的数字！！例如对于（894557.5 1025952.0）的第一行，我得到：（-82.355684 44.781735）。 - Rotail

负经度表示西经，正经度表示东经，而圣路易斯MO肯定在格林威治的西边 :) - Stef

你也可以通过谷歌进行验证。例如，数据框中的第一行（csv 中的第二行）在 csv 文件中具有地址 4124 Page Blvd。Google 地图为该地址提供了 -90.239668 / 38.650908 坐标，这与转换得到的 -90.239655 / 38.650896 相当接近，是合理的。 - Stef

我在想如何能够重现你的结果。我的计算机与你的一些设置似乎不同，这导致了结果上的差异。我尝试重置内核，但没有成功...我得到了相同的结果。 - Rotail