1、Cartopy vs Basemap
地球物理研究离不开画图，特别是画地图。Matplotlib-Basemap就是为了画题图而开发的，最初用来帮助和研究气候和天气预报的，但是长期以来它的API并不是特别好用，特别是在跨平台安装上。新一代路线图，准备用Cartopy替换basemap，未来basemap也将停止更新。
2020 年以后 Python 2.7 将停止更新，Basemap 会按照官方计划也迁移到 Cartopy 模块。所以，我们需要深入了解一下Cartopy和知道它有啥好。
首先， 官方文档说cartopy自带地理数据这一特征，大大简化了地图制图过程中准备数据的过程。
下面我们来试试这个库的效果，首先准备底图，Cartopy自带了一些基础数据，如：stock_img()是反映地形高程的背景。
另外，还有一些海岸线，湖泊边界之类的数据，直接add_feature函数可以控制 。
但有这些还不够，尤其是国界，最好用测绘局的，可以避免画错的问题，还有断层信息这些数据也需要自己给，add_geometries函数来搞定，测试了shp格式的本地矢量数据，添加进去没问题。
下面的图1就是小哥生成的一张地图，可以作为展示计算结果的底图。
图1 Cartopy的绘图效果
那么，添加信息怎么办呢？记住几个函数就行了plot和scatter是最常用的，让我们添加一些地震分布上去吧，scatter函数支持颜色和大小控制，添加地震后的效果如图2。这里大小与震级相关，颜色与震源深度相关。
图2 添加地震分布信息
有了这些可能还不够，还有有个图例，不着急，往下看，legend函数可以搞定图例，有时候自动生成的不好，还可以自定义，用*scatter.legend_elements设置即可，小哥试着加了一个图例，如图3所示。
图3 标注图例和南北地震带位置
上面几张图的代码，如下（为了看起来更加连贯，我就不一点点解释了，代码不复杂，很容易看懂，小哥也是从官网等渠道根据需要抄来的作业）：
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
from matplotlib.offsetbox import AnchoredText
datapath = Path(Path(catalog.__file__).parent,'data')
#1. 底图信息
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([72, 137, 10, 55])
ax.stock_img()
ax.coastlines()
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.RIVERS)
# 2.网格线
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
# 3.自定义信息
fname = Path(datapath, 'bou1_4l.shp')
f2name = Path(datapath, 'bou2_4l.shp')
faults = Path(datapath, 'gem_active_faults.shp')
ax.add_geometries(Reader(str(faults)).geometries(),
ccrs.PlateCarree(),facecolor = 'none',
edgecolor='red')
ax.add_geometries(Reader(str(f2name)).geometries(),
ccrs.PlateCarree(), facecolor = 'none',
edgecolor='gray', linestyle=':')
ax.add_geometries(Reader(str(fname)).geometries(),
ccrs.PlateCarree(), facecolor = 'none',
edgecolor='black')
# 4.地震分布
scatter = ax.scatter(data.longitude, data.latitude,
s= (0.2* 2 ** data.mag)**2,
c=data.depth / data.depth.max(), alpha=0.8,
transform=ccrs.PlateCarree())
# 5.标注图例
kw = dict(prop="colors", num= 6, fmt="{x:.0f} km",
func=lambda s: s*data.depth.max())
legend1 = ax.legend(*scatter.legend_elements(**kw),
loc="lower left", title="Depth")
ax.add_artist(legend1)
kw = dict(prop="sizes", num=5, color=scatter.cmap(0.7), fmt="M {x:.1f}",
func=lambda s: np.log2(np.sqrt(s)/0.2))
legend2 = ax.legend(*scatter.legend_elements(**kw),
loc="upper left", title="Mag")
# 6.标注范围
extent = (95, 108, 19.5, 44.5)
extent_box = sgeom.box(extent[0], extent[2], extent[1], extent[3])
ax.add_geometries([extent_box], ccrs.PlateCarree(), color='white',
alpha=0.5, edgecolor='black', linewidth=2)
# 7. 版权信息
SOURCE = 'GEOIST 2020'
LICENSE = 'MIT'
text = AnchoredText(r'$\mathcircled{{c}}$ {}; license: {}'
''.format(SOURCE, LICENSE),
loc=4, prop={'size': 12}, frameon=True)
ax.add_artist(text)
plt.show()
2、pygmt是GMT吗？
在地球物理学界，GMT的名气可谓无人不知。强大的矢量绘图能力和跨平台能力，为很多科研人员解决了高质量成图的困扰。
Python阵营这么红火，怎么少了GMT的接口呢？这不GMT也开始官方支持python了，即pygmt。
pygmt是GMT6.0之后的一个官方版本，以前当然也有很多类似的第三方接口，但是都不是官方的。
这里要强调的是pygmt还不是一个成熟、稳定的项目，还在开发中，我们测试一下效果即可。仿照上面的例子，先做底图后添加信息，绘图效果如图4所示，GMT那标志性的图框，熟悉吧，让我看看代码怎么写吧！
图4 pygmt的绘图效果
import pygmt as pg
fig = pg.Figure()
fig.basemap(region=region, projection=prj, frame=frame)
fig.coast(shorelines=sls, land="#666666", water="skyblue")
fig.plot(x=data.longitude,y=data.latitude,
sizes=0.01 * 2 ** data.mag,
color=data.depth / data.depth.max(),
cmap="viridis", style="cc", pen="black")
fig.savefig(filename)
代码不复杂，核心的就这几句话，注意一下上面的代码写法，是不是与matplotlib的pyplot接口用法很像，我想这就是pythonic GMT吧！但是，绘图过程会出现很多次dos的黑色弹出框，感觉不是太好。建议大家可以再等等，至少1.0版本出来再用。
一句话总结：今天写这个笔记，也是小哥自己的学习过程，现在Geoist里面混用了好几种地图绘制的库，后面会进一步根据最新的技术来封装Geoist自己的绘图函数，目标是尽量减少绘图任务和减小技术之间的耦合关系，Geoist开发也不想一下子提供太多的解决方案，而是期望做一套最好的即可，因为，坑趟过了，尽量不要让后人再掉进去。
网友评论