-
Notifications
You must be signed in to change notification settings - Fork 1
/
SC2SEI.py
281 lines (194 loc) · 11 KB
/
SC2SEI.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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu May 20 13:04:40 2021
@author: julio
based on: seiscomp_to_nordic.py
"""
# Programa para convertir un boletin de seiscomp acompañado del xml correspondie
# a formato nordic, este el formato utilizado por SEISAN en los S-Files
# los S-files constituyen la unidad de almacenamiento fundamental de SEISAN
# La definicion de nordic se puede hallar en el Apendice A del manual de SEISAN
#import os
from shutil import copyfile #biblioteca para manejo de archivo, especificamente copiar archivos
from seiscomp_xml_parser import xml_station_list #PARSER XML
import os # bliblioteca para utlidades de manejo de archivos
import sys
#funcion para agregar espacios en blanco a un string de una longitud determinada
#se agregan los espacios al principio del string s
def add_blank(s, length):
aux=s
if len(s)<length:
aux="{:>"+str(length)+ "}"
aux=aux.format(s)
return aux
#funcion para agregar espacios al final de un string de una longitud determinada
def add_blank_suffix(s, length):
aux=s
if len(s)<length:
aux="{:<"+str(length)+ "}"
aux=aux.format(s)
return aux
#funcion para determinar la posicion de un elemento dentro de una lista de listas
def find_in_sublists(lst, value):
for sub_i, sublist in enumerate(lst):
try:
return (sub_i, sublist.index(value))
except ValueError:
pass
raise ValueError('%s is not in lists' % value)
#sfile_name='/home/julio/Documents/project1/03-1925-42L.S202105'
#files_path='/home/julio/Documents/project1/'# directorio de trabajo
#seiscomp_file_name='insivumeh2021jfvi.txt'#nombre del boletin generado en seiscomp
#wav_file_name='insivumeh2021jfvi.mseed' #nombre del archivo wav al que apunta el s-file
#xml_file_name='insivumeh2021jfvi.xml' #xml generado por seiscomp, en este xml se tiene informacion de las arribos de ondas sismicas
def seiscomp_to_nordic(files_path,
seiscomp_file_name,
wav_file_name,
xml_file_name):
seiscomp_file=open(files_path+seiscomp_file_name,'r') #se lee el archivo del boletin
seiscomp_file_lines=seiscomp_file.readlines() # se leen las lineas del boletin de seiscomp
station_list_start=False # flag para determinar si se empiezan a leer picadas de fase
#station_list_header=False
station_list=[]# lista de estacione en el evento sismico
magnitude='Null'#valor temporal que indica que no se ha medido magnitud del evento sismico
depth=''
residual_RMS=''
for line in seiscomp_file_lines:# se itera a traves de las lineas del boletin con el fin de extraer informacion del evento
if line.find('Date')!=-1: #se extrae la fecha del evento sismico
year=line[27:31] # se lee el año
month=line[32:34]# se lee el mes
day=line[35:37]# se lee el dia
if line.find('Time')!=-1: # se lee la hora del evento sisimico
hour=line[27:29] # se lee la hora del evento sismico
minute=line[30:32] # se leen los minutos del evento sismico
seconds=line[33:37] # se leen los segundos del evento sismico
if line.find('Latitude')!=-1: # se lee la latitud del evento sismico
latitude=line[27:33]
if line.find('Longitude')!=-1: # se lee la longitud del evento sismico
longitude=line[27:33]
if line.find('Depth')!=-1: #se lee la profundidad del evento sismico
depth=line[31:33]
if line.find('Agency')!=-1:# se lee informacion de la agencia sismica
agency=line[27:30]
if line.find('Residual')!=-1: # se lee la informacion del rms
# print('test')
residual_RMS=line[29:33]
# print(line[29:33])
#if (line.find('ML')!=-1)and(line.find('MLv')==-1) and (magnitude=='Null'): # see lee la magnitud
# magnitude=line[14:18]
if ((line.find('MLv')!=-1)or (line.find('Mw')!=-1)) and (line.find('Preferred')==-1) and (magnitude=='Null'): # see lee la magnitud
magnitude=line[14:18]
if line.find('Public ID')!=-1: #id que sirve para identificar el origin utlizado por seiscomp para hallar el origen del evento sismico
public_ID=line[27:60]
station_list_header=False # se ha encontrado el header de las picadas de fase
if ((line.find('sta')!=-1)and(line.find('net')!=-1)and line.find('dist')!=-1 )and (len(station_list)==0):
station_list_start=True
station_list_header=True
if station_list_start and not station_list_header: # se extrae de las picadas de fase
station=line[4:9] #codigo de estacion
pick_hour=line[33:36] # hora de la llegada de fase
pick_minute=line[36:39] # minuto de la llegada de fase
pick_seconds=line[39:44]# segundo de la llegada de fase
res=line[45:49] #residual
pick_mode=line[50:51] #tipo de fase(P o S)
distance=line[16:19] #distacia del sismo
weight=line[53:56] #peso para la incerteza de la picada de fase
azimuth=line[20:23] #angulo azimutal
if (station not in station_list) and (station!=''):# no se tiene a la estacion en el listado de estaciones con picadas de fase
station_list.append([station, #0
pick_hour, #1
pick_minute, #2
pick_seconds, #3
res, #4
distance, #5
weight, #6
pick_mode, #7
azimuth]) #8
if station=='':
station_list_start=False
#################################################################################################################
# se utliza funcion creada en seiscomp_xml_parser para parsear el xml para extraer informacion de las estaciones con picada de fase
station_list_xml=xml_station_list(files_path+xml_file_name)
station_number=str(len(station_list_xml))# numero de estaciones encontradas
#depth=add_blank(depth,5)
station_number=add_blank(station_number,3)
residual_RMS=add_blank(residual_RMS, 4)
#magnitude=add_blank(magnitude)
#station_number=add_blank(station_number,3)
float_magnitude=float(magnitude)
#magnitude=f"{float_magnitude:.1f}"
magnitude="{:.1f}".format(float_magnitude)
magnitude=add_blank(magnitude,4)
#s_file=open(sfile_name,'r')
#s_file_lines=s_file.readlines()
residual_RMS_float=float(residual_RMS)
#residual_RMS=f"{residual_RMS_float:.1f}"
resdidual_RMS="{:.1f}".format(residual_RMS_float)
residual_RMS=add_blank(residual_RMS, 4)
###############################################################################################################
# se crea la primera linea del s-file, el header
month_aux=month # se arreglo bug para meses >=10
if month[0]=='0':
month_aux=' '+month[1]
s_file_header=' '+year+' '+month_aux+day+' '+hour+minute+' '+seconds+' L '+add_blank(latitude+'0',7)+add_blank(longitude+'0',8)+add_blank(depth+'.0',5)
s_file_header=s_file_header+' '+agency+station_number+residual_RMS+magnitude+'L'+agency
s_file_header=add_blank_suffix(s_file_header,79)+'1'
###############################################################################################################
# se crea la linea del s-file que contiene el nombre del wav-file(forma de onda)
seisan_wav_name=year+'-'+month+'-'+day+'T'+hour+':'+minute+':'+seconds[0:2]+'.MSEED'
#copyfile(files_path+wav_file_name,files_path+'/seisan_output/'+seisan_wav_name)
wav_file_location_line=' '+seisan_wav_name
wav_file_location_line=add_blank_suffix(wav_file_location_line,79)+'6'
sfile_name=day+'-'+hour+minute+'-'+seconds[0:2]+'L.S'+year+month
###############################################################################################################
# se crea el header que identifica al listado de picadas de fase sismicas
picks_header=' STAT SP IPHASW D HRMM SECON CODA AMPLIT PERI AZIMU VELO AIN AR TRES W DIS CAZ7'
###############################################################################################################
#se crea el archivo de salida es decir el s-file resultante, tambien se le agregan las lineas creadas arriba
#try:
# os.mkdir(files_path+'/seisan_output/')
#except:
# pass
copyfile(files_path+wav_file_name,files_path+seisan_wav_name)
s_file_output=open(files_path+ sfile_name,'w')
################################################################################################################
################################################################################################################
# inicia proceso de edicion de S-file
s_file_output.write(s_file_header+'\n')
s_file_output.write(wav_file_location_line+'\n')
s_file_output.write(picks_header+'\n')
#a=station_list_xml[0]
#se itera a traves de las picadas de fase para agregarlas al s-file
for station_pick in station_list_xml:
station=add_blank(station_pick[0],5)
channel=station_pick[1]
instrument=station_pick[2][0]
hour=station_pick[6]
minute=station_pick[7]
seconds=station_pick[8]
phase=station_pick[9]
quality_indicator='I'# revisar Apendice A del manual de SEISAN pag 541..
station_pick_line=' '+add_blank_suffix(station_pick[0],5)+instrument+channel+' '+quality_indicator+phase
station_pick_line=station_pick_line+' '+hour+minute+seconds+'0 '
station_pick_line=add_blank_suffix(station_pick_line,80)#+azimuth[:6]+' '
#station_pick_line=add_blank_suffix(station_pick_line,63)+time_residual+' '+distance+' '
s_file_output.write(station_pick_line+'\n')
s_file_output.write(add_blank(' ',80)+'\n')
s_file_output.close()
print(sfile_name+' '+xml_file_name)
#return s_file, seisan_wav_name
#root_path='./SC2SEI/'
root_path=sys.argv[1]
seicomp_file_list=os.listdir(root_path)
i=0
for seiscomp_file in seicomp_file_list:
i+=1
#print(i)
#print(seiscomp_file)
if seiscomp_file.find('xml')!=-1:
seiscomp_file_name=seiscomp_file.split('.')[0]
seiscomp_to_nordic(root_path,
seiscomp_file_name+'.txt',
seiscomp_file_name+'.mseed',
seiscomp_file)