內(nèi)容摘要：書(shū)接上文，繼續(xù)聊Geoist中的catalog模塊的相關(guān)功能，給大家介紹如何獲取EQ目錄數(shù)據(jù)和進(jìn)行相應(yīng)的分析與質(zhì)量控制QC。主要是給EQ目錄分析感興趣和給剛?cè)胄械男“讉兛磭D！
都兩口氣了還沒(méi)寫(xiě)完午磁，因?yàn)榧纫貓D，又要敲字毡们，實(shí)在不容易，再接再厲昧辽，我們繼續(xù)下面4個(gè)部分內(nèi)容衙熔。
7、重復(fù)事件檢測(cè)
What搅荞？一個(gè)EQ目錄里面如果記錄了兩個(gè)相同的事件红氯，就叫重復(fù)事件。
Why咕痛？啥怎么還能記錄重了痢甘！因?yàn)橛袝r(shí)候EQ定位不準(zhǔn)，茉贡。
How塞栅？原理知道了，下面我們看怎么整：
真心說(shuō)腔丧，USGS目錄還不錯(cuò)放椰，也可能是5級(jí)EQ起步，重復(fù)的可能性不大愉粤，圖1按照兩個(gè)EQ相差16s砾医，距離在100km之內(nèi)的事件，列為疑似重復(fù)記錄衣厘。跑了一篇檢測(cè)結(jié)果才3對(duì)如蚜。當(dāng)然，根據(jù)不同的數(shù)據(jù)標(biāo)準(zhǔn)還可以改呦！
8错邦、目錄截取
What探赫？截取就是取一部分。
Why兴猩？因?yàn)橛袝r(shí)候分析目錄需要一部分就夠了期吓，數(shù)據(jù)多了計(jì)算慢，所以咱的截一塊倾芝。
How讨勤？前面講了一種用Filter的方法，原理我們都知道了晨另，下面我們看看還能怎么整：
代碼如下：
from geoist.catalog import Selection as sel
p = [(90.,20.),(90.,40.),(105.,40.),(105.,20.),(90.,20.)]
db3 = sel.AreaSelect(db2,p)
上面用到了AreaSelect函數(shù)潭千，方法是先定義一個(gè)多邊形p，然后借尿，將目錄db2中的數(shù)據(jù)選出來(lái)刨晴，賦值給一個(gè)新的db3的目錄對(duì)象。
9路翻、目錄平滑
What狈癞？目錄平滑是指通過(guò)一定算法，將每個(gè)空間范圍內(nèi)都給一個(gè)發(fā)震概率或其它量茂契，有時(shí)候一些EQ風(fēng)險(xiǎn)評(píng)估模型蝶桶，需要空間的連續(xù)函數(shù)來(lái)作為輸入。
Why掉冶？EQ怎么還要平滑呢真竖！不行咱們盤(pán)幾下？厌小！ 反正無(wú)論怎么樣恢共，我就想要來(lái)空間上的連續(xù)函數(shù)。
How璧亚？我們這里的平滑采用高斯函數(shù)讨韭，原理大家都知道，從上面選擇的db3開(kāi)始涨岁，下面我們看怎么整：
from geoist.catalog import Smoothing as sm
from geoist.catalog import CatUtils as ct
from geoist.catalog import Exploration as exp
x1,y1,z1 = exp.GetHypocenter(db3)
P = ct.Polygon()
P.Load(p)
wkt = ct.XYToWkt(P.x, P.y)
xsm, ysm, asm = sm.SmoothMFD(db3, 1., wkt, Delta=0.5)
cfg1 = {'Bounds': [90., 20., 105., 40.],
'FigSize': [10., 12.],
'Background': ['none',[0.9,0.8,0.6],[0.5,0.8,1.]],
'Grid': [5., 5.]}
m1 = mapt.GeoMap(cfg1)
m1.BasePlot()
m1.MeshPlot(xsm, ysm, asm)
#m1.AreaPlot(P.x, P.y, Set=['y',0.5,'k',1])
#m1.PointPlot(xsm, ysm, Set=['o','b',2,1], Label='Grid')
m1.PointPlot(x1, y1, Set=['o','g',5,1], Label='全部')
m1.DrawGrid()
m1.Title('川滇地區(qū)EQ目錄高斯平滑')
m1.Show()
上面代碼運(yùn)行沒(méi)問(wèn)題拐袜，會(huì)輸出下圖：
圖中綠色小點(diǎn)，是這個(gè)區(qū)域內(nèi)的EQ震中梢薪，基本哪里EQ多何荚，哪里顏色紅无午，但是這些背景顏色是空間上連續(xù)函數(shù)，這有時(shí)候計(jì)算模型是十分需要的格式荠医。
10、去余震（Decluster）
What？大EQ之后，會(huì)誘發(fā)一系列余震，把這些挑出來(lái)就是去余震過(guò)程恒序。
Why？EQ事件可以分成前震谁撼、主震歧胁、余震。一般主震是那個(gè)震級(jí)最大的厉碟，因?yàn)橛嗾鸬陌l(fā)生主要?dú)w因于主震的誘發(fā)喊巍，也就是說(shuō)這類(lèi)EQ活動(dòng)不應(yīng)該算到正常的EQ活動(dòng)里面去。所以箍鼓，在有些分析EQ活動(dòng)的場(chǎng)景崭参，需要去掉余震事件。那怎么去呢款咖？因?yàn)橹髡饌€(gè)頭大何暮，余震除了震級(jí)小，就是距離主震震中距離近铐殃，時(shí)間短海洼，所以，辦法就是富腊，確定一個(gè)規(guī)則來(lái)篩選贰军。
How？去余震也叫decluster過(guò)程蟹肘，主要還是按照EQ事件和震級(jí)來(lái)算，震級(jí)越大我們要考慮更長(zhǎng)時(shí)間范圍內(nèi)俯树，主震周邊的EQ為余震帘腹。原理知道了，下面我們看怎么整：
這里面用了一個(gè)叫GardnerKnopoff的去余震方法许饿，代碼如下：
from geoist.catalog import Declusterer as declus
dbm, log1 = declus.WindowSearch(dbusgs, WinFun= declus.GardnerKnopoff, WinScale=1)
dbm.Info()
為去余震的目錄信息：
Number of Events: 5383
Year Rage: (1970,2019)
Magnitude Rage: (5.0,7.9)
Latitude Rage: (15.017,54.979)
Longitude Rage: (70.02,134.95)
Depth Rage: (0.0,588.0)
去掉可能是余震的事件阳欲，信息如下：
Number of Events: 2714
Year Rage: (1970,2019)
Magnitude Rage: (5.0,7.9)
Latitude Rage: (15.017,54.979)
Longitude Rage: (70.02,134.95)
Depth Rage: (0.0,588.0)
5383到2714是不是少了不少啊陋率！
11球化、震級(jí)轉(zhuǎn)換
What？一個(gè)EQ目錄里面如果用了不同的震級(jí)標(biāo)度瓦糟，這時(shí)候要進(jìn)行一些統(tǒng)計(jì)就比較麻煩筒愚，因?yàn)闃?biāo)準(zhǔn)不統(tǒng)一。
Why菩浙？不同的EQ巢掺，以震中與記錄臺(tái)站之間的距離句伶，可以分為遠(yuǎn)震和近震等，這時(shí)候?yàn)榱朔奖懵降恚赡艹霈F(xiàn)不同的震級(jí)標(biāo)度考余，這就會(huì)引起統(tǒng)一個(gè)EQ，兩家單位定出來(lái)的EQ震級(jí)不一樣轧苫，類(lèi)似于英里和公里楚堤。
How？嚴(yán)格來(lái)講EQ震級(jí)之間不能轉(zhuǎn)換含懊，但是有時(shí)候?yàn)榱擞?jì)算和分析也必須要轉(zhuǎn)身冬，下面我們看怎么整：
度量一個(gè)EQ的震級(jí)方法有有不少，所以震級(jí)單位也有很多绢要，比如：mb吏恭，ml，ms重罪，mw樱哼。嚴(yán)格來(lái)講這些震級(jí)之間是不能轉(zhuǎn)換的，但是如果你不是很較真剿配，方法還是有的搅幅。這類(lèi)方法叫經(jīng)驗(yàn)公式法，就是大概率下是差不多的呼胚，有限應(yīng)用目的沒(méi)問(wèn)題茄唐。
下面我們將mb震級(jí)轉(zhuǎn)為mw，代碼如下：
from geoist.catalog import Selection as sel
from geoist.catalog import MagRules as mr
sel.MagConvert(dbm,'*',['mb','Mb'],'Mw',mr.mb_Mw_Lin_DiGiacomo2015)
dbm.Info()
輸出信息：
Converting ['mb', 'Mb'] to Mw: 1236 events found
對(duì)比轉(zhuǎn)換前后蝇更，兩個(gè)同樣的事件信息沪编，dbm.Events[99]：
轉(zhuǎn)換前：
{'Id': 'usp00003gf',
'Log': 'earthquake',
'Location': [{'Year': 1973,
'Month': 9,
'Day': 21,
'Hour': 13,
'Minute': 48,
'Second': 32.2,
'Latitude': 18.834,
'Longitude': 120.65,
'Depth': 20.0,
'SecError': None,
'LatError': None,
'LonError': None,
'DepError': None,
'LocCode': None,
'Prime': False}],
'Magnitude': [{'MagSize': 5.2,
'MagError': None,
'MagType': 'mb',
'MagCode': None}]}
轉(zhuǎn)換后：
{'Id': 'usp00003gf',
'Log': 'earthquakeMAGCONV(None:mb);',
'Location': [{'Year': 1973,
'Month': 9,
'Day': 21,
'Hour': 13,
'Minute': 48,
'Second': 32.2,
'Latitude': 18.834,
'Longitude': 120.65,
'Depth': 20.0,
'SecError': None,
'LatError': None,
'LonError': None,
'DepError': None,
'LocCode': None,
'Prime': False}],
'Magnitude': [{'MagSize': 5.39,
'MagError': 0.0,
'MagType': 'Mw',
'MagCode': None}]}
這里將mb=5.2轉(zhuǎn)換完事mw=5.39.
Mw震級(jí)物理意義明確，直接與能量對(duì)應(yīng)年扩，能統(tǒng)一成Mw震級(jí)是最好的蚁廓，小哥的Geoist庫(kù)里面為大家準(zhǔn)備了一個(gè)中國(guó)大陸及周邊的數(shù)據(jù)。文件名為mwcat1900utc.csv厨幻，我們可以叫cnmw相嵌，詳情請(qǐng)參考資料3。如果要使用這個(gè)格式的目錄導(dǎo)入庫(kù)况脆，請(qǐng)參考下面代碼：
catname = 'cnmw'
localpath = qc.pathname+'\\mwcat1900utc.csv'
db6 = cat.Database(catname)
header = ['Year', 'Month','Day','Hour','Minute','Second','Latitude','Longitude','MagSize','Depth','Log']
db6.Import0(localpath, Header = header, Delimiter= ',', flag = False)
不知道大家是否發(fā)現(xiàn)饭宾，這個(gè)header定義很靈活，只要這部分信息沒(méi)錯(cuò)格了，任何格式目錄都是可以輕松玩轉(zhuǎn)的看铆。
一句話(huà)總結(jié)：哇，EQ目錄分析好有意思啊盛末，大家反映上次寫(xiě)的太科普性湿，小哥這不發(fā)奮學(xué)習(xí)了好幾天纬傲，發(fā)揚(yáng)動(dòng)手能力強(qiáng)的優(yōu)點(diǎn)，咱也能攢個(gè)分析工具出來(lái)肤频！現(xiàn)在工具也有了叹括，數(shù)據(jù)也有了，屏幕前奮發(fā)圖強(qiáng)的小盆友們宵荒，趕快動(dòng)手吧汁雷，下一次開(kāi)組會(huì)，PPT里面的圖可以用用geoist的catalog工具包氨取侠讯！
這次講完了:)
參考資料：
1、USGS目錄暑刃，https://earthquake.usgs.gov/fdsnws/event/1/
2厢漩、CENC目錄，中國(guó)EQ臺(tái)網(wǎng)發(fā)布的標(biāo)準(zhǔn)目錄岩臣，小哥做了內(nèi)網(wǎng)地址溜嗜，http://10.2.14.222/catalog
3、中國(guó)矩震級(jí)目錄架谎，詳情參考文章：Cheng, J., Y. Rong, H. Magistrale, G. Chen, and X. Xu (2017). An Mw‐based historical earthquake catalog, Bull. Seismol. Soc. Am. 107 2490-2500.
4炸宵、GEOIST工具包，你關(guān)心的信息全在這里https://github.com/igp-gravity/geoist
