番外篇4_Geoist之地图的画法

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就是小哥生成的一张地图，可以作为展示计算结果的底图。

上面几张图的代码，如下（为了看起来更加连贯，我就不一点点解释了，代码不复杂，很容易看懂，小哥也是从官网等渠道根据需要抄来的作业）：

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那标志性的图框，熟悉吧，让我看看代码怎么写吧！
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(),