在地球科学领域，需要画很多地形相关的图件，平时使用比较多的是ArcGIS和Adobe系列软件，还有免费的GMT。在读博期间我也修过用GMT 画图的课，用GMT 画图的确很好。但是由于平时用的不多，每次画图我都要一点点查找命令，一般很熟的人几分钟搞定的底图，我要折腾半天。虽然最终可以搞定，但是效率实在不高。虽然Julia 已经渐渐流行了，但是用了Python 十年的我还是比较喜欢Python。
众所周知，Python 很强大，一般用来下载数据，处理数据，并最终成图，用Python画的图一般会很美观。想着如果可以用Python 画图那多有趣，有了底图，我再将Python 处理出来的数据画在上面，就会得心应手。
有了探索的心，就弱弱地尝试过几次，当时用的Python 2.7 版本，最后也没有成功。近来，我更新了系统，安装了Python3.6. 一个偶然的机会，成功地用Python 画了全球地形图。与大家在此分享，如果想探索更多，请自行修改代码。
测试环境
操作系统：Ubuntu 18.04Python 版本：Python 3.6.9
安装 basemap 模块
这个网上可以搜到很多种方法，我尝试过很多次，最终发现这个可以成功：
conda install basemap
代码实现
因为conda 安装的原因，需要在用时改变一下环境变量，在代码前面加一段：
import os
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
from mpl_toolkits.basemap import Basemap
完整代码样例
本例子基于文献(1)
#!/usr/bin/env python
#!/usr/bin/env python
import matplotlib.pyplot as plt
import csv
from io import StringIO
import numpy as np
import math
import requests
from csv import DictReader
################ For the environment value change for basemap
import os
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
from mpl_toolkits.basemap import Basemap
################ read data from the downloaded files
key='4.5_month.csv'
### read eq data
aa=np.genfromtxt(key,skip_header=0,names=True,dtype=None,delimiter=',',invalid_raise = False,deletechars="~!@#$%^&*()-=+~\|]}[{';: /?.>,<.", case_sensitive=True,usecols=(0,1,2,3,4))
### read original data
time=[]
lats=[]
lngs=[]
depth=[]
mags=[]
mags_adjust=[]
for i in range(0,len(aa)):
time.append(str(aa[i][0]))
lats.append(float(aa[i][1]))
lngs.append(float(aa[i][2]))
depth.append(float(aa[i][3]))
mags.append(float(aa[i][4]))
mags_adjust.append(2.2**float(aa[i][4]))
'''
################ Alternatively, we can download earthqauke data from USGS website directly instead of the above section
#URL for USGS for Earthquake data
DATA_URL = 'https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv'
print("Downloading", DATA_URL)
resp = requests.get(DATA_URL)
print(resp)
quakes = list(DictReader(resp.text.splitlines()))
lngs = [float(q['longitude']) for q in quakes]
lats = [float(q['latitude']) for q in quakes]
mags = [float(q['mag']) for q in quakes]
mags_adjust = [2.2**float(q['mag']) for q in quakes]
'''
################ plot basemap of earth!!! it is very nice!!!
plt.figure(figsize=(15, 9))
earth = Basemap()
earth.bluemarble(alpha=0.42)
earth.drawcoastlines(color='#555566', linewidth=1)
plot1=plt.scatter(lngs, lats, mags_adjust, c='red',alpha=0.5, zorder=10, edgecolors='none')
marker1=plt.scatter([],[],s=2.2**(min(mags)),c='red',alpha=0.5,)
marker2=plt.scatter([],[],s=2.2**(min(mags)+1.),c='red',alpha=0.5,)
marker3=plt.scatter([],[],s=2.2**(min(mags)+2.),c='red',alpha=0.5,)
labels=[str(min(mags)),str(min(mags)+1),str(min(mags)+2)]
plt.legend([marker1,marker2,marker3],labels,ncol=4,labelspacing=1.2, handlelength=3, frameon=True, borderpad=1,title='Earthquake size',loc='lower left', markerscale=1.,scatterpoints=1,fontsize=10)
plt.title("M4.5+ earthquakes in recent 30 days")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.xticks(range(-180,180,10))
plt.yticks(range(-90,90,10))
plt.grid()
plt.savefig('usgs-4.5quakes-recent_30_days.png', dpi=350)
效果图
网友评论