pyproj 2：CRS与高斯坐标转换

CRS

CRS即坐标参照系(coordinate reference system)，是pyproj中的核心类，可以对坐标进行更复杂的描述，并且对EPSG提供了更多的构造方式，以4326代表的WGS84为例，方法如下

# 下面函数都是生成WGS84坐标系统
from pyproj import CRS
crs = CRS.from_epsg(4326)
crs = CRS.from_string("epsg:4326")
crs = CRS.from_proj4("+proj=latlon")
crs = CRS.from_user_input(4326)
1
2
3
4
5
6

一般来说，我国的基本比例尺地形图较多采用高斯克吕格投影，根据地图比例尺不同，又采用不同度数的分带。

高斯投影

其中，高斯六度带和UTM分带的计算方法相似，但起始经度不同，其东半球的分带号可表示为

$$N=\lfloor L/6\rfloor +1$$

三度带则为

$$N=\lfloor (L+1.5)/3\rfloor$$

一般来说，比例尺大于1:1万的地图，均采用高斯三度带投影。

北京的经纬度是$(116.425,39.906)$，则其在高斯三度带投影中为39号，在投影坐标系统中的名称为 Xian_1980_3_Degree_GK_Zone_39，对应的EPSG号为2363，则其转换方法为

from pyproj import Transformer

# 这里采用CRS类，也可以直接用字符串来表示
crs = CRS.from_epsg(2363)
tf = Transformer.from_crs("epsg:4326", crs) 
 
lat = 39.906
lon = 116.425
x, y = tf.transform(lat, lon)
print("x:", x, "y:", y)
# x: 4419252.202604582 y: 39450831.09905656
1
2
3
4
5
6
7
8
9
10
11

其输出值中，经度所对应的$y$值明显大得离谱，已知纬度1°对应111km，经度由于要乘上一个三角函数，所以必然小于此值，而输出的$y$值却达到了39450km这个级别，这都快赶上地球周长了。

但稍加思索就会看到，这个39正好是带号，故而最终得到的平面直角坐标为$(4419252.202604582,450831.09905656)$，这样看来就合理了很。

proj的转换函数tf.transfrom不仅支持单个经纬度的转换，也可以输入numpy数组，

lon = np.linspace(116, 117, 10)
lat = np.linspace(30, 40, 10)
x, y = tf.transform(lat, lon)
1
2
3

在西安80坐标系中，其高斯三度带所对应的带号从25到45，相应地EPSG编号从2349到2369，总计大约垮了60多度的经线范围。

高斯六度带所对应的带号从13到23，相应地EPSG编号从2327到2337。