美文网首页程序员
使用Python 画地形图与地震分布!湖南才是最安全的!

使用Python 画地形图与地震分布!湖南才是最安全的!

作者: 59aa689f8c96 | 来源:发表于2019-09-27 14:34 被阅读0次

在地球科学领域,需要画很多地形相关的图件,平时使用比较多的是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)

效果图 

相关文章

网友评论

    本文标题:使用Python 画地形图与地震分布!湖南才是最安全的!

    本文链接:https://www.haomeiwen.com/subject/csfauctx.html