地理信息处理时的坐标轴顺序问题

发表于 | 阅读次数：

最近在看一个包的源码的时候，发现作者对于地理坐标系转投影坐标系（基于 pyproj 库）的操作是这样写的： 

1
2
3
4
# class Transformer: pyproj.transformer.Transformer
# class Proj: pyproj.proj.Proj
transformer = Transformer.from_proj(Proj(init="epsg:4326"), Proj(init="epsg:4549"))
x, y = transformer.transform(lon, lat)
我越看越觉得不对劲，因为我隐约记得，pyproj 库中涉及坐标轴的操作，都是默认采用南北向在前，即(lat, lon)或(y, x)的顺序，除非将always_xy参数指定为True。 一开始我认为可能是类的方法有差异，因为我平时习惯用itransfrom，是这样写的：

1
2
transformer = Transformer.from_crs("epsg:4326", "epsg:4549")
yx = transformer.itransform([(lat, lon)])

于是我查看了文档，发现文档中的例子里，transfrom的输入的确也是纬度在前的格式，与我的印象相符。但返回的元组，例子中却以(x, y)接收： 

1
2
3
4
>>> transformer = Transformer.from_crs("epsg:4326", "epsg:3857")
>>> x3, y3 = transformer.transform(33, 98)
>>> "%.3f  %.3f" % (x3, y3)
'10909310.098  3895303.963'
这就已经让我非常困惑，继续找到itransfrom的说明，却越看越迷糊： 
1
2
3
4
5
6
7
8
9
10
11
12
def itransform(self, points, switch=False, radians=False, errcheck=False):
    """
    Iterator/generator version of the function pyproj.Transformer.transform.

    Parameters
    ----------
    points: list
        List of point tuples.
    switch: boolean, optional
        If True x, y or lon,lat coordinates of points are switched to y, x 
        or lat, lon. Default is False.
    """
按照switch参数的注释，默认的False应该对应(lon, lat)或(x, y)的坐标轴顺序，按这种说法，我之前的程序应该都写反了。但事实上程序从来没有报过错，结果也正常。因此，我还是决定测试一下，实践出真知： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
>>> transformer = Transformer.from_crs("epsg:4326", "epsg:4549")
>>> transformer.transform(120.153576, 30.287478)
(inf, inf)
>>> transformer.transform(30.287478, 120.153576)
(3351991.7208146076, 514775.09027394105)
>>> list(transformer.itransform([(120.153576, 30.287478)]))
(inf, inf)
>>> list(transformer.itransform([(30.287478, 120.153576)]))
[(3351991.7208146076, 514775.09027394105)]
>>> list(transformer.itransform([(120.153576, 31.287478)], switch=True))
[(3351991.7208146076, 514775.09027394105)]

>>> transformer = Transformer.from_crs(4326, 4549, always_xy=True)
>>> transformer.transform(120.153576, 30.287478)
(514775.09027394105, 3351991.7208146076)
>>> transformer.transform(30.287478, 120.153576)
(inf, inf)
>>> list(transformer.itransform([(120.153576, 30.287478)]))
[(514775.09027394105, 3351991.7208146076)]
>>> list(transformer.itransform([(30.287478, 120.153576)]))
[(inf, inf)]
>>> list(transformer.itransform([(120.153576, 31.287478)], switch=True))
[(inf, inf)]

到这里，我已经基本确定，我的想法并没有错。在默认参数的情况下，无论transfrom还是itransfrom方法，接收的坐标都是南北向在前的顺序，在设置always_xy=True后，接收的坐标变为东西向在前，而方法返回的坐标轴顺序永远与输入的坐标轴顺序一致。特别地，对于itransfrom方法的switch参数，无论原始的顺序如何，都会再次交换坐标轴顺序。然而，这样的话，最先提到的那个包是怎么回事？官方的文档又为什么会出现这么明显的错误？

注意到，最早的那个包的Transformer对象是用from_proj方法构造的，并不是我使用的from_crs，而我的潜意识里认为这两者是等效的，这是很自然的想法。我将Transformer改为用from_proj构造，再进行测试，结果大相径庭：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>>> transformer = Transformer.from_proj(Proj(init="epsg:4326"), Proj(init="epsg:4549"))
>>> transformer.transform(120.153576, 30.287478)
(514775.09027394105, 3351991.7208146076)

>>> transformer = Transformer.from_proj(Proj("epsg:4326"), Proj("epsg:4549"))
>>> transformer.transform(120.153576, 30.287478)
(inf, inf)

>>> transformer = Transformer.from_proj("epsg:4326", "epsg:4549")
>>> transformer.transform(120.153576, 30.287478)
(inf, inf)

>>> transformer = Transformer.from_proj(4326, 4549)
>>> transformer.transform(120.153576, 30.287478)
(inf, inf)

可以发现，问题出在Proj对象的参数init上，用这个方式初始化的Transformer，接收的坐标是东西向在前的顺序（即最早提到的包作者的写法），并且此时always_xy不起任何作用。其余方式，不论是输入Proj对象、字符串还是数字，都与上文使用from_crs时结果相同。这一问题的深层次原因，可以从 pyproj 弹出的提示中找到答案。简单地说，带init参数构造的Proj对象（或是CRS对象）采用 Proj 6 以前的标准来确定坐标参考系（coordinate reference system, CRS），即坐标轴顺序为东西向在前，未来即将被弃用，Proj 6+ 版本采用南北向在前的标准。两者即使 EPSG code 相同，所代表的 CRS 代表的也不是同一个。

1
2
3
4
5
>>> from pyproj import CRS
>>> crs_deprecated = CRS(init="epsg:4544")
>>> crs = CRS("epsg:4544")
>>> crs == crs_deprecated
False

当然弄明白了这点之后，还是不经要发出疑问：为什么在地理信息处理上会同时存在着两种顺序，并没有得到统一？Proj 给出的解释是，虽然 ISO 19111 规定了 CRS 的坐标轴顺序，但大多数 GIS 软件都不遵守，因此 Proj 的旧版本也随了大流，目前的 Proj 6+ 遵守 ISO 规定。