1嘹屯、Cartopy vs Basemap
地球物理研究離不開畫圖抚垄，特別是畫地圖。Matplotlib-Basemap就是為了畫題圖而開發(fā)的呆馁，最初用來幫助和研究氣候和天氣預(yù)報(bào)的，但是長(zhǎng)期以來它的API并不是特別好用阴挣，特別是在跨平臺(tái)安裝上纺腊。新一代路線圖，準(zhǔn)備用Cartopy替換basemap誓沸，未來basemap也將停止更新壹粟。
2020 年以后 Python 2.7 將停止更新，Basemap 會(huì)按照官方計(jì)劃也遷移到 Cartopy 模塊洪添。所以雀费，我們需要深入了解一下Cartopy和知道它有啥好。
首先盏袄， 官方文檔說cartopy自帶地理數(shù)據(jù)這一特征，大大簡(jiǎn)化了地圖制圖過程中準(zhǔn)備數(shù)據(jù)的過程炭菌。
下面我們來試試這個(gè)庫的效果逛漫，首先準(zhǔn)備底圖，Cartopy自帶了一些基礎(chǔ)數(shù)據(jù)克握，如：stock_img()是反映地形高程的背景枷踏。
另外，還有一些海岸線停团，湖泊邊界之類的數(shù)據(jù)，直接add_feature函數(shù)可以控制 佑稠。
但有這些還不夠，尤其是國界捆蜀，最好用測(cè)繪局的幔嫂，可以避免畫錯(cuò)的問題，還有斷層信息這些數(shù)據(jù)也需要自己給锰茉，add_geometries函數(shù)來搞定切心，測(cè)試了shp格式的本地矢量數(shù)據(jù)，添加進(jìn)去沒問題。
下面的圖1就是小哥生成的一張地圖定鸟，可以作為展示計(jì)算結(jié)果的底圖。
那么啼县，添加信息怎么辦呢沸久？記住幾個(gè)函數(shù)就行了plot和scatter是最常用的，讓我們添加一些地震分布上去吧子刮，scatter函數(shù)支持顏色和大小控制窑睁，添加地震后的效果如圖2。這里大小與震級(jí)相關(guān)橱赠，顏色與震源深度相關(guān)箫津。
有了這些可能還不夠宰啦，還有有個(gè)圖例赡模，不著急惕耕，往下看，legend函數(shù)可以搞定圖例司澎，有時(shí)候自動(dòng)生成的不好，還可以自定義谚殊，用*scatter.legend_elements設(shè)置即可蛤铜，小哥試著加了一個(gè)圖例，如圖3所示剿干。
上面幾張圖的代碼置尔，如下（為了看起來更加連貫，我就不一點(diǎn)點(diǎn)解釋了榜轿，代碼不復(fù)雜朵锣，很容易看懂，小哥也是從官網(wǎng)等渠道根據(jù)需要抄來的作業(yè)）：
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.網(wǎng)格線
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.標(biāo)注圖例
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.標(biāo)注范圍
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. 版權(quán)信息
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嗎助析？
在地球物理學(xué)界椅您，GMT的名氣可謂無人不知。強(qiáng)大的矢量繪圖能力和跨平臺(tái)能力掀泳，為很多科研人員解決了高質(zhì)量成圖的困擾西轩。
Python陣營這么紅火藕畔，怎么少了GMT的接口呢？這不GMT也開始官方支持python了庄拇，即pygmt。
pygmt是GMT6.0之后的一個(gè)官方版本溶弟，以前當(dāng)然也有很多類似的第三方接口瞭郑，但是都不是官方的。
這里要強(qiáng)調(diào)的是pygmt還不是一個(gè)成熟擒权、穩(wěn)定的項(xiàng)目阁谆，還在開發(fā)中，我們測(cè)試一下效果即可剖效。仿照上面的例子裳凸，先做底圖后添加信息劝贸，繪圖效果如圖4所示，GMT那標(biāo)志性的圖框梦湘，熟悉吧件甥，讓我看看代碼怎么寫吧！
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)
代碼不復(fù)雜瓣颅，核心的就這幾句話譬正，注意一下上面的代碼寫法檬姥，是不是與matplotlib的pyplot接口用法很像健民，我想這就是pythonic GMT吧贫贝！但是，繪圖過程會(huì)出現(xiàn)很多次dos的黑色彈出框稚晚，感覺不是太好。建議大家可以再等等筑辨，至少1.0版本出來再用幸逆。
一句話總結(jié)：今天寫這個(gè)筆記，也是小哥自己的學(xué)習(xí)過程楚昭，現(xiàn)在Geoist里面混用了好幾種地圖繪制的庫拍顷，后面會(huì)進(jìn)一步根據(jù)最新的技術(shù)來封裝Geoist自己的繪圖函數(shù)，目標(biāo)是盡量減少繪圖任務(wù)和減小技術(shù)之間的耦合關(guān)系尿贫，Geoist開發(fā)也不想一下子提供太多的解決方案，而是期望做一套最好的即可庆亡，因?yàn)槔谈澹犹诉^了，盡量不要讓后人再掉進(jìn)去彰亥。
