Python画任意两点间的地形剖面图

灭火器

之前一直不知道怎么画任意两点间的地形剖面，这两天试了一下，记录在这里。
问题1：地形数据从哪里来？
如果是处理 WRF 数据，那么 wrfout 中就有地形数据。除此之外，可以自己到网络上下载地形起伏数据。下面这篇博客简要介绍了常用的数据：
全球地形起伏数据总结：https://blog.seisman.info/global-relief-models/
这里使用的是 ETOPO2v2c 的全球地形数据，格式为 nc，网格为等经纬度网格。下载地址为：
https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO2/ETOPO2v2-2006/

问题2：怎么获取任意两点间连线的地形？
给出一个点，我们可以利用等距网格的数据插值得到这个点的高度的数值。那么给出任意两点间连线上一系列的点，对每个点进行这样的操作，就能获取所有点（也就是整条线）上的地形高度。这里利用 scipy.interpolate 模块的 interpn 函数进行插值，这个函数专门用来处理规则网格的数据-->散点的插值。或者也可以自己写双线性插值的函数，结果与 interpn 是一致的。

代码和效果图如下。其中地形的水平分布填色图摘自帖子：
Python绘制全球地形和中国地形：http://bbs.06climate.com/forum.php?mod=viewthread&tid=95503&highlight=%B5%D8%D0%CE
这里直接用的 numpy 数组进行的插值。如果你会 xarray 的话，可以考虑使用 metpy 库中的 cross_section 函数，一行代码即可出结果，并且两点间的连线使用的是测地线。参考连接如下：
https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.cross_section.html#metpy.interpolate.cross_section

更新
之后在坛友对方的凤飞飞的帮助下，实现了极地投影的地形剖面图。代码进行了如下修改：
（1）使用 ETOPO2v2g 作为数据源，因为其经度是首尾相连的。如果使用 v2c 的话，则画不出正常的 contour 图。
（2）极地处为了保证连线的正常，使用pyproj库生成的测地线。
（2）修改了经纬度的标签。
效果图如下：

两份代码贴在附件中。最后吐槽一下，用 cartopy 画极地真的坑，标签都要手动加，感觉还是 NCL 方便……

对方的凤飞飞

高级,真的学习,真的高级

qqqql

感谢分享。学习学习

likailing2008

厉害厉害，这个可能是地理学的大佬

对方的凤飞飞

请问楼主,处理极地投影这种投影似乎不可行呀,有办法解决吗

灭火器

对方的凤飞飞 发表于 2020-9-27 11:53
请问楼主,处理极地投影这种投影似乎不可行呀,有办法解决吗

试了下，用cartopy的话，对于南极的投影，我连水平的地形图都画不出来，时间又长，contour level还是错的。我看你帖子说十几秒就能画出来，我真傻眼了，不如你教教我吧。

对方的凤飞飞

灭火器 发表于 2020-9-27 15:27
试了下，用cartopy的话，对于南极的投影，我连水平的地形图都画不出来，时间又长，contour leve ...

可以呀,我好像已经解决了,我可以私发给你

灭火器

对方的凤飞飞 发表于 2020-9-27 15:58
可以呀,我好像已经解决了,我可以私发给你

试了半天，发现是用 ETOPOv2c 画不出，换成 v2g 就好了。

对方的凤飞飞

灭火器 发表于 2020-9-28 09:58
试了半天，发现是用 ETOPOv2c 画不出，换成 v2g 就好了。

那是我换了数据,怎么样还凑合吧

yrovl

太厉害了{:eb502:}{:eb502:}