本月第一轮讲授航图实践课程。虽说名称叫实践，但以往的执行以纸笔绘图练习为主。完成了障碍物A型图、精密进近地形图之后，隐约感觉大家上路之后进度可能很快，于是在航路图这节打算上点强度：认认真真地把经纬度投影出来做航路图的练习。对于非测绘或者地理科班的学生而言，地图投影可能一直以来都是听个概念，很少动手实际体验这个过程——即使是地图学专业，我记忆里似乎也少有手工展绘坐标点的经历，只要在软件或者代码里设置好转出转入的坐标系也就够用了。
问题产生：投影坐标轴/顺序定义
这节课的练习里要绘制一块2°×2°的航路图局部，首先面临的一个问题就是图廓线的绘制。无论选择什么投影，这个球面上的整齐区域对应的平面图廓一定不应该是个正方形——虽然对于没有制图基础的人来说正方形是最符合直觉的。这时候，就需要根据制图区域的四至范围计算图廓四个角点的投影坐标。甚至可以更粗略地，计算左下角点和右上角点的坐标来确定一个矩形框作为图廓。
方便易懂的epsg.io网站是我手动转换坐标时的首选工具。利用该网站，可以把WGS 84 (EPSG:4326)的坐标转换为投影坐标。这里值得考虑的投影有Lambert、UTM、Gauss-Krüger1。兰伯特因为不同地区的中央经线和标准纬线不同，在众多搜索结果中不好找合适的，就没有在课堂上使用。这里先把转换工具支持的ESRI:102012 (WGS 1984 Lambert Asia)结果在表格最后列出供对比。UTM和GK的选择只需要看经度范围，简单易行，此处列出111°中央经线六度带的投影WGS 84 / UTM 49N (EPSG:32649) CGCS2000 / GK 19 (EPSG:4497)。
| 坐标系 | 左下 22N 110E | 右上 24N 112E | 左上 | 右下 |
| EPSG:4497 | 19396734, 2434138 | 19601755, 2655649 | 19398244, 2655649 | 19603265, 2434138 |
| EPSG:32649 | 396775, 2433164 | 601714, 2654587 | 398285, 2654587 | 603224, 2433164 |
| ESRI:102012 | 539685, 2819337 | 734587, 3065085 | 525045, 3049050 | 755070, 2835819 |
有了平面坐标，就可以分别在横纵两个坐标轴上分别作差得到制图区域的实地跨度，进而按照比例尺转化为图面距离即图廓的尺寸了。其实至此理论上已经没有问题了，但这个很稳的事情再次由于课堂条件的限制出现了问题：时间有限、难度不能太高，课堂上我只算了GK19的前两个(也就是左下和右上)角点，引出了一个关键问题，哪个是长、哪个是宽？或者说哪个是北、哪个是东？
问题发展：两个反直觉的事实
复盘来看上文表中完整的3个投影下各4个角点共12个坐标，这根本不是问题。最直观的，左上和右上的Y坐标相等，所以Y坐标是上下/南北方向。深一层的，GK投影的坐标如果不算带号，东西数字位数少、南北数字位数多，算上两位带号反过来2，示例点在6度的19带，19开头的长数字是带号，也就是X是东西方向。因此，制图区域东西长205km、南北长221km，为了放在A4纸21cm宽的幅面上还能留有一定边距，比例尺取1:120万。
课堂条件限制下，只算了两个点所以最直观的思路失效，而第一次选用6度GK分带方案则导致了不熟悉19是个带号，进而使专业背景带来的经验规律失效。题外话，最近武测地信教材以丛书的形式在测绘出版社出版，我采购了地信教程第2版和投影原理与实践两种，摘录其中一个基本事实备忘：我国6°投影带13~23带、3°投影带24~45带。
其实专业经验也并没有完全失效，我意识到了测绘习惯以北为X方向，然而这点残存的经验规律配合epsg.io网站的错误验证导致了错误的产生。
Coordinate system: Cartesian 2D CS. Axes: northing, easting (X,Y). Orientations: north, east. UoM: m.
EPSG:4497 CGCS2000 / Gauss-Kruger zone 19
Source: OGP3, as cited at https://epsg.io/4497
如上，epsg.io据称是原样照搬的EPSG数据库显示EPSG:4497是个平面直角坐标系，X是北、Y是东、单位是米。这让我当时有点确信地认为南北短、东西长，而实际上网站提供的计算结果(上表)是东在前、北在后。这是第一个反直觉的事实，即X是东且带号为19。
第二个反直觉的事实是不专业的，课堂上同学(当然不仅仅是学生，普通大众突然被问到这个问题通常也会这么想)提供的支持错误结论的证据，大概是大众地理知识加上不怎么深入的思考带来的错误推论。这个直觉结论的表面层次是，平时见到的世界地图都是东西长、南北短，因为经度±180°、加起来一共360°，而纬度只有±90°、加起来一共180°，360＞180，所以我们的制图区域也是东西长、南北短。这个直觉的错误在于，世界地图的形状和本次制图的形状没有关系，不能类比。同时，如果说这个平时见到的投影是朴素的Plate Carrée (简易圆柱)投影，那相同经纬度跨度的区域就应该是正方形。
这个错觉的深层次是，地球是一个梨形形状，赤道半径大于极地半径，所以相同角度的弧长东西长、南北短。这个错觉是因为忽略了赤道附近之外的地区，东西方向扫过弧长对应的半径是在赤道半径的基础上乘以纬度的余弦，这个余弦带来的差异远比地球长轴短轴带来的差异明显：6356/6378=0.996, arccos0.996=4.76°。也就是说，在纬度高于5°的地区，余弦的效应更明显，即东西向的单位弧长短于南北向。
问题纠正：展绘点位与坐标刻度的不一致
由于以上条件限制下的错觉，在课堂上发生了轻微翻车的一幕。有两位同学表示，根据epsg.io转换的地物投影坐标折合比例尺后展绘的点位，与第一步(错觉下错误判断的)图廓线刻度不一致。但这比起竟无一人是男儿的沉默大多数还是令人欣慰的。其中一位更是坚持要理清，否则如果后续出于练习目的，将按刻度绘制的要素与按投影坐标和比例尺绘制的要素共冶一炉，可能导致空间关系的错乱。虽然担心的错乱完全可以通过只采用刻度描点或者只采用投影换算描点(应该不会，毕竟手工太麻烦了)得以避免，但这个思维之清晰和一丝不苟还是难能可贵的。
其实上一节弧长的逻辑就是发现错误后同学给出反对错误结论论据的思路。从测绘专业的角度以及我本人的习惯，发现上述问题后我的第一反应是，坐下来把所有的点坐标全部列出，画出草图检查。而学员直接按照常识的角度搜索一度对应的实际距离，得知一纬度对应实地111公里，而一经度对应实地111cosφ公里，比如纬度40°的一经度对应实地85公里。所以投影转换工具给出的坐标是东西在前、南北在后。
讨论：软件工具中的投影坐标定义与实现
事已至此，epsg.io网站上提供的EPSG:4326转EPSG:4497肯定是存在某些问题的，因为上文引用的文字说明和表格计算结果是自相矛盾的——CGCS的GK 19带定义是X表示北，但计算结果中X却表示东。于是，忍不住去看了开源的epsg.io网站的实现。
这是一个Python Flask应用，同时也有ChromeOS的版本放在同一个代码库中。从javascript/src/transformpage.js中可以看到页面监听的一系列变化事件都指向了epsg.io.TransformPage.prototype.transform_函数，而这个函数在校验完输入后用goog.net.Jsonp向一个指向'/trans'节点的服务请求结果数据。再转到app.py中的该服务节点实现，根据解析出的坐标系用osgeo.osr模块创建了CoordinateTransformation对象，采用该转换器类的TransformPoint方法完成了坐标转换并写入到响应结果中。与此同时，关于坐标系的信息存储在一个sqlite数据库中，在app.py中的多个与坐标系编号相关节点中均可见用sqlite3模块创建的数据库连接和游标查询读取结果的语句。
也就是说，epsg.io这个坐标转换工具网站声称的从EPSG数据库及其他可靠来源拉取坐标系资料是真实的，EPSG对CGCS的GK投影的描述(X北、Y东)也是合理的。但是这个工具网站在底层实际调用了OSGeo中的OSR也就是GDAL (geospatial data abstract library)中的OGR Spatial Reference的坐标转换实现(的Python接口，原函数库是C++)。这倒也正常，毕竟GDAL也足够可靠和权威，但问题是即使是这样的库，对坐标顺序的定义问题的解决也是近几年的事情。在这个文档中提到，GDAL从3.0版本(19年5月)开始才支持按坐标系定义的顺序返回结果，在这之前，该库默认以所谓的传统GIS顺序(先经后纬、先东后北)返回结果，无论坐标系如何定义。所以这次被坑的原因就在于选用的坐标转换工具网站调用的是旧版本的底层代码，无视了坐标系定义的顺序，导致与网页从EPSG数据库拉取的描述信息不一致。
当然，太多年不看代码，并不是一开始就找到上述思路和结论的。之前写Python代码的时候常用的地图投影库是pyproj，其底层似乎是Proj.44。所以这次遇到状况第一反应是之前用这个库时收到过一个升级重大变化提醒。于是打开Jupyter测了一下，坐标顺序完全遵守原始定义，才让我下定决心去翻工具网站的实现一探究竟。题外话，现在AI写代码是真方便，连帮助文档都不用查，在百度搜“pyproj投影例子”直接就生成带注释的代码，让我这个面向Stack Overflow编程流派感到老了。
>>> import pyproj
>>> wgs84 = pyproj.CRS('epsg:4326')
>>> gk19 = pyproj.CRS('epsg:4497')
>>> tfm = pyproj.Transformer.from_crs(wgs84, gk19)
>>> tfm.transform(22, 110)
(2434138.0838125898, 19396734.053191964)
>>> utm49 = pyproj.CRS('epsg:32649')
>>> tfm2 = pyproj.Transformer.from_crs(wgs84, utm49)
>>> tfm2.transform(22, 110)
(396775.35957092687, 2433164.428653589)讨论：中等尺度区域的投影选择
回到课程和练习的问题本身。从航行的角度来说，航路图用UTM/GK与现行规章是有一定偏差的，这只是课堂和练习限制下的妥协。如前文所述，采用Lambert投影会增加理解难度，表中的数据也表明制图区域不在中央经线上会导致图幅上下框线与纸面边缘不平行、不水平；同时，纬线的弧线半径很大也难以采用作图工具绘制，其纸笔操作难度也大幅上升；再者，自定义的标准纬线和中央经线的兰伯特投影转换缺乏像epsg.io网站这样方便的工具，虽然用代码不难实现，但对于非制图、非计算机专业学生而言也是令人窒息的挑战。以上三点与课程让学生熟悉图面要素、要素符号、关联的AIP资料表格查询与解读等目标相背离。
民用航空仪表航路图及区域图采用正轴等角双标准纬线圆锥投影。航路图标准纬线24°和40°，中央经线106°；区域图根据实际图幅范围确定标准纬线和中央经线。
WM-TM-2021-002, §4.2.3
采用正轴等角双标准纬线圆锥投影。
MH/T 4019-2012, §8.2 §9.2
（航路图）应使用直线接近大圆的正形投影（例如，兰伯特正形投影）。
（区域图）应采用直线接近大圆的正形投影。兰伯特圆锥正形投影非常适合作为绘制此类图的依据。投影类型不需要在图上注明。
ICAO Doc 8697, §7.7 §7.8
建议：应采用直线接近大圆的正形投影。
ICAO Annex 4, §7.4.1 §8.4.1
从制图的角度，即使是航路图的局部，比例尺也较小。如本次练习中2°的范围在A4纸上比例尺就超过了百万分一。虽然制图中从来都没有明确的大小比例尺的界限，但小于1:100万的比例尺着实不能算大。前述投影教材的GK章节中提到50万和100万比例尺的图中可以绘制经纬网，似乎又在暗示GK投影也可以有如此小比例尺的应用。ESRI的指导手册里认为，兰伯特正形圆锥投影使用的范围是大陆/大洋、地区/海域、大比例尺，高斯-克吕格和横轴墨卡托则是后两项，差别只是少了大陆/大洋一项。网络上的一些讨论也认为如果制图区域和一个UTM网格的大小相近，则对UTM投影的使用是合适的。总之，GK/UTM绘制2°/200km尺度的范围肯定不是最佳，但也绝不是完全不合适，属于无可无不可的状态。
所以本次练习的决定以及轻微出错只能说是一种妥协的产物，正路在于在繁忙的工作任务中把航图实践课程的制图部分升级到数字化和自动化的解决方案。