- 2021-12-23 22:38:16
前些日子的一篇博文中, 在谈到分析环状分子的puckering坐标时, 简单提及了用gnuplot绘制球面投影的方法. 那里所用的几种投影的变换公式都比较简单, 不涉及复杂的求解. 实际中广泛应用的一种投影是Mollweide投影, 我们见到的世界地图大多使用这种投影方法绘制. 这种投影的变换公式复杂一些, 要对每个经纬度求解一个非线性方程才能得到变换后的平面坐标. 求解时可以使用牛顿迭代法, 具体参考维基百科中的说明.
在gnuplot中实现牛顿迭代法并不困难, 但如果想将其封装成一个函数, 那就必须绕一些弯弯. 这是因为gnuplot所支持的函数只能使用一条语句, 而不能直接在其中使用判断, 循环等复合语句. 而我们要用的牛顿迭代法, 那是必须要进行迭代的. 怎么解决这个矛盾呢? 我想到的方法有两种, 一种基于eval, 先将要执行的语句拼接成字符串, 然后用eval执行拼好的字符串. 这种方法并不能算是完美地实现了函数, 因为没有办法直接获得返回值, 而且必须使用全局变量, 但可以算是比较通用的一种解决方法, 因为任何语句都可以这么做. 另一种方法就是将迭代改成写递归, gnuplot只支持尾递归, 递归深度也有限制(好像是255), 所以不是任何迭代都可以这么做, 好在我们需要的牛顿迭代方法是可以改写成尾递归形式的.
实现了Mollweide投影方法之后, 我们就可以根据球面上每个点的经纬度将其转换为平面坐标, 如果再将每个地点的一些性质也一起绘制出来, 就成为地图了. gis(地理信息系统)程序做的就是类似的事, 研究地球物理的也经常做这种事, GMT是这个领域比较常用的程序. 用gnuplot做这种事不专业, 但可以学到gnuplot的一些高级用法, 虽然免不了编写脚本.
在编写gnuplot脚本的时候, 我才发现它竟然支持使用utf-8字符作变量, 所以使用中文字符也就没有问题了. 有了这个”大”发现, 我一度想将脚本中的所有变量都用中文定义, 从而让脚本看起来像中文语言的样子, 但试了试之后就放弃了, 觉得这并不是好的做法. 用中文写公式虽然可行, 但已经废弃很久了. 至于原因, 读读晚清李善兰翻译的全中文微积分课本, 你就理解了.
我不知道现在中学的解题格式是什么样子的了, 我中学时候的理科考试, 试卷是中文的, 但答题时也没有人完全用中文写公式, 大多是先说明一下符号的含义, 后面就直接用那些符号了. 所用符号大致是英文字母和希腊字母, 下标什么的可以用中文, 如令时间为t, 初速度为v0, 则其位移为s=v0*t之类, 和教科书上的写法类似. 你真要完全写成中文, 我觉得还更难记, 而且也耗费更多时间. 所以呢, 一种编程语言支持使用任意字符做变量是好的, 但这不意味着我们要用它做完全的中文编程. 真要那么做的话, 建议你学学文言编程, 格调还更高一些. 所以, 我们可以在脚本中有限制合理地使用一些中文字符, 希腊字符, 如λmin=起始经度=-180; λmax=终止经度=180; λ0=中心经度=0; dλ=经度步长=10之类, 起到注释作用, 让脚本变得更易读易懂是可以的, 但完全的中文编程, 既没有必要, 也浪费时间.
还是让我们回到地图的绘制吧.
就让我们绘制全球的海拔高度和海洋深度图吧, 这是我们熟悉的. 为此, 我们需要全球的高程和海深数据. 公开的这类数据不是很多, 比较全的有两份, 一份来自GEBCO, 一份来自NASA. 前者2021版数据的分辨率为15弧秒, 涵盖了海拔和水深, 数据完整. 后者是汇编多个数据集而成的, 分辨率为1弧秒, 只有海拔, 大致涵盖了全球80%的陆地区域. 这类全球数据都很大的, 动辄上G, 多使用二进制格式, 处理使用起来都比较麻烦, 需要进行转换.
绘制海拔高度还需要合适的颜色映射方案. 常用颜色方案有GMT自带的那些, M_Map修改过的, 还有一个我觉得比较漂亮的 , 它们的特点是看起来直观, 与实际情况比较符合.
如果需要添加大陆的边界线, 可以使用GMT或Natural Earth的海岸线数据, 后者有些地方并不精确. 如果需要世界各国及其行政区划的数据, GMT带有DCW的数据, 也可以试试gadm的数据, 前者最多到省级边界, 后者包含了直到县区的边界. 但要注意, 这两份数据给出的国界线未必符合各国的领土主张, 中国的尤甚, 只给出了大陆部分的边界, 没有包含港澳台, 南海部分, 藏南部分也不合适. 真要绘制涉及中国的部分, 建议使用GMT中文社区整理的中国地理空间数据集.
对于更小的区域, GEBCO数据的分辨率就不够了, 可以使用NASADEM的数据, 虽然有些地方不全.
上面大致把我了解到的信息概述了一下, 因为我不是专门研究gis的, 所以有些地方没弄对也可以理解. 这种东西, 用gnuplot绘制不是不可以, 但GMT可能更方便, 毕竟它是专门用于绘制地图的. 当然, 也有一些matlab的工具, 如M_Map: A mapping package for Matlab.
最后来看看我绘制的一些结果吧. 先是全球的高程海深数据, 几种颜色映射方案看起来区别并不大. 更好的效果要使用Hillshade方法进行晕染, 可惜我还没完成.
再看看中国部分
南中国海特写
最后看看我家乡吧, 因为我更熟悉那里.
