-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot_events_from_2010_Z-800_two_stations.py
113 lines (95 loc) · 4.88 KB
/
plot_events_from_2010_Z-800_two_stations.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
"""Maciej Barnaś, [email protected]
Program plots stations and events from 2010 year, with local magnitude greater than 1.2, below the 600 m b.s. l in
Bobrek coal mine (data from IS-EPOS platform).
Last edit: 2017-10-02"""
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.basemap import pyproj
from datetime import datetime
from def_import_Bobrek import import_events, import_stations
from numpy import arange
print(__doc__)
Bobr = import_events(file='BOBREK_catalog_01.2010_ML1.2_Z-600.mat') # Importing events from Matlab file
stations = import_stations()
PL2000VI = pyproj.Proj("+init=EPSG:2177") # ETRS89 / Poland CS2000 zone 6
# Convert geographical coordinates of the stations to Poland CS2000 zone 6 coordinates
for lab, row in stations.iterrows():
stations.loc[lab, 'Y_2000'], stations.loc[lab, 'X_2000'] = PL2000VI(row['Longitude'], row['Latitude'])
stations_two = stations[stations['Code'].str.contains('S013|S015')]
# Convert geographical coordinates of the events to Poland CS2000 zone 6 coordinates
for lab, row in Bobr.iterrows():
Bobr.loc[lab, 'Y_2000'], Bobr.loc[lab, 'X_2000'] = PL2000VI(row['Long'], row['Lat'])
# Adding column with event number
for lab, row in Bobr.iterrows():
Bobr.loc[lab, 'nr'] = Bobr.loc[lab, 'ID'][17:]
Bobr[['nr']] = Bobr[['nr']].apply(pd.to_numeric)
events_2010 = Bobr['nr'] >= 932 # I checked on the IS-EPOS platform, that this event is first
# after 2010 Jan 19 06:09:48.0 (date of installing some sensors), it would be difficult, to
# check it in Python, because of 'Time' column - it doesn't contain dates, but big floats,
# so converting it to dates takes a lot of time
Bobr1 = Bobr[events_2010] # Events from 2010 Jan 19 06:09:48.0 to the end of the data
events_2010_below_800 = Bobr1['Z'] < -800
Bobr2 = Bobr1[events_2010_below_800]
events_2010_below_800_ML12 = Bobr2['ML'] >= 1.2
Bobr3 = Bobr2[events_2010_below_800_ML12]
fig = plt.figure()
ax1 = fig.add_subplot(121, projection='3d')
plt.tight_layout()
ax1.scatter(Bobr3['Y_2000'], Bobr3['X_2000'], Bobr3['Z'], s=30*Bobr3['ML'], c=Bobr3['ML'], cmap='rainbow')
ax1.scatter(stations_two['Y_2000'], stations_two['X_2000'], stations_two['Elevation'], marker='^', s=40, color='r')
for lab, row in stations_two.iterrows():
ax1.text(stations_two.loc[lab, 'Y_2000']+3, stations_two.loc[lab, 'X_2000']+3, stations_two.loc[lab, 'Elevation']+3,
stations_two.loc[lab, 'Code'])
# ax1.set_xlabel('Latitide [\u00b0]')
# ax1.set_ylabel('Longitude [\u00b0]')
ax1.set_xlabel('Y [m]')
ax1.set_ylabel('X [m]')
ax1.set_zlabel('Z [m]', labelpad=20)
ax1.tick_params(axis='z', pad=12)
# max_lat = 50.3716
# min_lat = 50.3432
# max_long = 18.8982
# min_long = 18.8329
min_X_2000 = 5.579e6 + 250
max_X_2000 = 5.579e6 + 1750
min_Y_2000 = 6.56e6 + 1000
max_Y_2000 = 6.56e6 + 2500
x_ticks = arange(5.579e6 + 500, 5.579e6 + 2000, 500)
y_ticks = arange(6.56e6 + 1000, 6.56e6 + 3000, 500)
ax1.set_xlim([min_Y_2000, max_Y_2000])
ax1.set_ylim([min_X_2000, max_X_2000])
ax1.set_xticks(y_ticks)
ax1.set_yticks(x_ticks)
ax2 = fig.add_subplot(122)
#D = ax2.scatter(Bobr2['Y_2000'], Bobr2['X_2000'], c=Bobr2['ML'], cmap='rainbow', s=30)
#ax2.scatter(stations_two['Y_2000'], stations_two['X_2000'], marker='^', s=40, color='r')
#D = ax2.scatter(Bobr3['Long'], Bobr3['Lat'], c=Bobr3['ML'], cmap='rainbow', s=30)
D = ax2.scatter(Bobr3['Y_2000'], Bobr3['X_2000'], c=Bobr3['ML'], cmap='rainbow', s=30)
#ax2.scatter(stations_two['Longitude'], stations_two['Latitude'], marker='^', s=40, color='r')
ax2.scatter(stations_two['Y_2000'], stations_two['X_2000'], marker='^', s=40, color='r')
#plt.plot([18.8670, 18.8717], [50.3547, 50.3546], c='y', linewidth=5)
#plt.plot([18.8670, 18.8717], [50.3515, 50.3515], c='y', linewidth=5)
for lab, row in stations_two.iterrows():
ax2.text(stations_two.loc[lab, 'Y_2000']+3, stations_two.loc[lab, 'X_2000']+3, stations_two.loc[lab, 'Code'])
# for lab, row in stations_two.iterrows():
# ax2.text(stations_two.loc[lab, 'Longitude'], stations_two.loc[lab, 'Latitude'], stations_two.loc[lab, 'Code'])
ax2.set_xlabel('Y [m]')
ax2.set_ylabel('X [m]')
# ax2.set_xlabel('Longitude [\u00b0]')
# ax2.set_ylabel('Latitude [\u00b0]')
cbar = fig.colorbar(D, orientation='horizontal')
cbar.set_label('Local magnitude')
# ax2.set_xlim([min_long, max_long])
# ax2.set_ylim([min_lat, max_lat])
ax2.set_xlim([min_Y_2000, max_Y_2000])
ax2.set_ylim([min_X_2000, max_X_2000])
ax2.set_xticks(y_ticks)
ax2.set_yticks(x_ticks)
ax2.grid()
plt.figtext(0.5, 0.03, 'Figure 3. Seismic events with magnitude higher than 1.2 and location below 800 m b.s.l. '
'(points) - color responds magnitude; two stations chosen for further calculations - red'
'triangles with labels; on the left - 3D plot, on the right - 2D map',
horizontalalignment='center', fontsize=11, style='italic', backgroundcolor='w')
print(stations_two)
plt.show()