-
Notifications
You must be signed in to change notification settings - Fork 0
/
3_NND.py
131 lines (104 loc) · 5.38 KB
/
3_NND.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#!python2.7
'''
Created on April 10th, 2019
- compute nearest-neighbor distance (NND) between all event pairs (see equ. 1 in Zaliapin & Ben-Zion, 2013)
- test -plot histogram of NNDs: Figure 4c in Zaliapin & Ben-Zion, 2013
output: 'data/%s_NND_Mc_%.1f.mat'%(dPar['catName'], dPar['Mc'])
which is a python dictionary with:
{ 'aNND' : aNND, - nearest neighbor space-time magnitude distance
'aEqID_p' : np.array - ID of the parent event
'aEqID_c' : np.array - ID of the child event
}
### dNND only contains clusterred events
TODO:
- constrain Mc, b and D independently through statistical analysis of the actual data
@author: tgoebel
'''
#------------------------------------------------------------------------------
import matplotlib as mpl
mpl.use('Agg') # turn off interactive plot
import matplotlib.pyplot as plt
import numpy as np
import scipy.io
import os
#------------------------------my modules--------------------------------------
import src.clustering as clustering
from src.EqCat import *
eqCat = EqCat( )
#print dir( dataUtils)
#=================================1==============================================
# dir, file, params
#================================================================================
dir_in = '%s/data'%( os.path.expanduser('.'))
file_in= 'hs_1981_2018_all.mat'
#file_b = '%s_b_Mc_D.txt'%(fileIn.split('.')[0])
dPar = { 'aMc' : np.array( [3.1,3.2,3.3,3.4,3.6,3.7,3.8,3.9]),#np.array([3.0, 4.0]),np.array( [2.0, 2.5, 3.0, 3.5]
# fractal dimension and b for eq. (1)
'D' : 1.6, #1.6 TODO: - these values should be contrained independently
'b' : 1.0,
#=================plotting==============
'eta_binsize' : .3,
'xmin' : -13, 'xmax' : 0,
}
#=================================2==============================================
# load data, select events
#================================================================================
eqCat.loadMatBin( os.path.join( dir_in, file_in))
print( 'total no. of events', eqCat.size())
eqCat.selectEvents( dPar['aMc'][0], None, 'Mag')
#eqCat.selector( tmin, tmax, 'Time')
print( 'no. of events after initial selection', eqCat.size())
#=================================1==============================================
# to cartesian coordinates
#================================================================================
# two ways to do the distance comp: 1 project into equal distance azimuthal , comp Cartersian distance in 3D
# 2 get surface distance from lon, lat (haversine), use pythagoras to include depth
eqCat.toCart_coordinates( projection = 'aeqd')
for dPar['Mc'] in dPar['aMc']:
print( '-------------- current Mc:', dPar['Mc'], '---------------------')
# select magnitude range
eqCat.selectEvents( dPar['Mc'], None, 'Mag')
print( 'catalog size after MAG selection', eqCat.size())
# this dictionary is used in module: clustering
dConst = {'Mc' : dPar['Mc'],
'b' : dPar['b'],
'D' : dPar['D']}
#==================================2=============================================
# compute space-time-magnitude distance, histogram
#================================================================================
dCluster = clustering.NND_eta( eqCat, dConst, correct_co_located = True)
###histogram
aBins = np.arange( -13, 1, dPar['eta_binsize'])
aHist, aBins = np.histogram( np.log10( dCluster['aNND'][dCluster['aNND']>0]), aBins)
aHist, aBins = np.array(list(zip(*zip(aHist, aBins))))# cut to same length
# correct for binsize
aHist /= dPar['eta_binsize']
# to pdf (prob. density)
aHist /= eqCat.size()
#=================================3==============================================
# save results
#================================================================================
import scipy.io
NND_file = 'data/%s_NND_Mc_%.1f_HD.mat'%( file_in.split('.')[0], dPar['Mc'])
# NND_file = 'data/%s_NND_Mc_%.1f.mat'%( file_in.split('.')[0], dPar['Mc'])
print( 'save file', NND_file)
scipy.io.savemat( NND_file, dCluster, do_compression = True)
#=================================4==============================================
# plot histogram
#================================================================================
fig, ax = plt.subplots()
#ax.plot( vBin, vHist, 'ko')
ax.bar( aBins, aHist, width =.8*dPar['eta_binsize'], align = 'edge', color = '.5', label = 'Mc = %.1f'%( dPar['Mc']))
ax.plot( [-5, -5], ax.get_ylim(), 'w-', lw = 2, label = '$N_\mathrm{tot}$=%i'%( eqCat.size()))
ax.plot( [-5, -5], ax.get_ylim(), 'k--', lw = 2, label = '$N_\mathrm{cl}$=%i'%( dCluster['aNND'][dCluster['aNND']<1e-5].shape[0]))
ax.legend( loc = 'upper left')
ax.set_xlabel( 'NND, log$_{10} \eta$')
ax.set_ylabel( 'Number of Events')
ax.grid( 'on')
ax.set_xlim( dPar['xmin'], dPar['xmax'])
plt.show()
plotFile = 'plots/%s_NND_hist_Mc_%.1f.png'%( file_in.split('.')[0], dPar['Mc'])
print( 'save plot', plotFile)
plt.savefig( plotFile)
plt.show()
plt.clf()