Python地理坐标系和投影坐标系转换

0 相关名词

在开始之前,有必要了解一下相关名词:

  1. 地心地固坐标系(Earth-Centered, Earth-Fixed,ECEF),简称地心坐标系。
  2. 地理坐标系统(Geographic Coordinate System,GCS)1,坐标系是地心坐标系,用经纬度表示球面上的点。
  3. 世界大地测量系统(World Geodetic System, WGS),比如WGS84,是一种地理坐标系统,用于全球定位系统(GPS)。
  4. 投影坐标系统(Projection Coordinate System,PCS)2,,在二维平面上用表示位置。
  5. 通用横轴墨卡托投影(Universal Transverse Mercator,UTM),是一种投影方法。

关于地理坐标系和投影坐标系更详细的解释可以查看这篇文章:你必须知道的地理坐标系和投影坐标系

地理坐标系统有不同的基准和方法,比如:Xian_1980,Beijing_1954,WGS_1984等。投影坐标系统也有不同的基准和方法,UTM和UPS等。每一个地理坐标系统(GCS)和投影坐标系统(PCS)都有一个独特的EPSG代码,代码可在 EPSG 网站查询。

有一篇介绍Pyproj进行地理投影坐标系转换的文章3,但不够全面。其中提到arcgis网站上查询 地理坐标系投影坐标系 的方法很实用但不全。

1 地理和投影坐标系统相互转换

整理使用Python的第三方库 Pypro4 转换经纬度表示的地理坐标系统和(或千米等其他单位)表示的投影坐标系统。

Pypro模块共有两个函数:

函数 描述
test() 运行模块测试
transform(p1, p2, x, y, z=None, radians=False) 用法:x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False),将在坐标系统p1下的点(x1, y1, z1)转换到p2坐标系统下

1.1 使用EPSG Code转换

转换经纬度到的投影坐标系统;转换一个投影坐标系统到另一个投影坐标系统;反向转换,把投影坐标系统上的点转换到地理坐标系统:

>>> p1 = pyproj.Proj(init='epsg:26915') # 一个投影坐标系统EPSG Code>>> p2 = pyproj.Proj(init='epsg:26715') # 另一个投影坐标系统EPSG Code>>> x1, y1 = p1(-92.199881,38.56694) # 投影到EPSG Code为26915的投影坐标系统>>>> '%9.3f %11.3f' % (x1,y1)'569704.566 4269024.671'>>> x2, y2 = pyproj.transform(p1,p2,x1,y1) # 转换一个投影坐标系统到另一个投影坐标系统>>> '%9.3f %11.3f' % (x2,y2)'569722.342 4268814.027'>>> '%8.3f %5.3f' % p2(x2,y2,inverse=True) # 反向转换' -92.200 38.567'

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

用元组传入多个点

>>> lats = (38.83,39.32,38.75) # 所有纬度组成的元组>>> lons = (-92.22,-94.72,-90.37) # 所有精度组成的元组>>> x1, y1 = p1(lons,lats) # 转换经纬度到投影坐标系统>>> x2, y2 = pyproj.transform(p1,p2,x1,y1) # 转换一个投影坐标系统到另一个投影坐标系统>>> lons, lats = p2(x2,y2,inverse=True) # 反向转换
  • 1
  • 2
  • 3
  • 4
  • 5
  • 1
  • 2
  • 3
  • 4
  • 5

1.2 使用基准名称转换

除了使用EPSG Code之外,还可以显示指定坐标系统名称

>>> p1 = pyproj.Proj(proj='latlong',datum='WGS84') # WGS84,GPS使用的地理坐标系统,EPSG Code为4326>>> x1 = -111.5; y1 = 45.25919444444>>> p2 = pyproj.Proj(proj='utm',zone=10,datum='NAD27') # 投影坐标系统NAD27 / UTM zone 10N,EPSG Code为26710>>> x2, y2 = pyproj.transform(p1, p2, x1, y1)>>> '%s %s' % (str(x2)[:9],str(y2)[:9])'1402285.9 5076292.4'

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

2 总结

如果仅仅是需要把经纬度转换为米为单位的话,可直接按照下边来:

>>> p1 = pyproj.Proj(init='epsg:4326') # 定义数据地理坐标系 WGS84>>> p2 = pyproj.Proj(init='epsg:3857') # 定义转换投影坐标系>>> x, y = pyproj.transform(p1, p2, lon, lat) # lon 和lat 可以是元组
  • 1
  • 2
  • 3
  • 1
  • 2
  • 3
(0)

相关推荐