Example # 01: Comparing two Hercules incidents

Prepared by: Naeem Khoshnevis (nkhshnvs@memphis.edu, khoshnevis.naeem@gmail.com)

Disclaimer: These results are only presented for illustrating different features of tsprocess. The runs are some random simulations. They should not be used for any other purposes.

See Setting Up Environment section to set up the environment.

In [1]:
# import project module
from tsprocess import project as pr
In [2]:
# Let's start a new project
p1 = pr.Project("project_1")

Each project will connect to a database with project name + "_db.sqlite" database (in this example: project_1_db.sqlite). If the database is exist, it will use that database, otherwise new database will be generated.

In [3]:
# Let's take a look at a summary of the project database.
p1.database_summary()

Database name: project_1_db.sqlite
Database size: 0.012288 MB.
Number of items: 1
List of database trackers: ['project_1_dbtracker']

* database tracker is used to track records of different incidents inside the database. You can ignore it.

In [4]:
# Each project needs a source hypocenter. These tests runs are conducted for 2014 Lahabra, CA, earthquake.
p1.add_source_hypocenter(33.933, -117.916, 5.1)

Please note that these folders are prepared according to Incidents preparation instructions.

Here are the description.txt of these folders.

# folder herclules_1 description.txt:    
incident_name             = hercules101    
incident_type             = hercules     
inputfiles_parameters     = inputfiles/parameters.in     
source_hypocenter         = 33.933, -117.916, 5.1     
hr_comp_orientation_1     = 39.9      
hr_comp_orientation_2     = 129.9      
ver_comp_orientation      = down      
incident_unit             = m
# folder herclules_2 description.txt:         
incident_name             = hercules102         
incident_type             = hercules        
inputfiles_parameters     = inputfiles/parameters.in        
source_hypocenter         = 33.933, -117.916, 5.1       
hr_comp_orientation_1     = 39.9        
hr_comp_orientation_2     = 129.9      
ver_comp_orientation      = down      
incident_unit             = m
In [5]:
# let's add incidents to the project
In [6]:
h_folder_1 = 'hercules_1'
h_folder_2 = 'hercules_2'
In [7]:
p1.add_incident(h_folder_1)
p1.add_incident(h_folder_2)

Adding incident to the project connects the incident folder to project and loads its stations location. However, it does not load any record. Inside the description.txt file there is a incident_name variable. It should be a unique value. Any attemp to load new incident with used names will be denied. For example, if we want to add incident for h_folder_1 again, it will be denied.

In [8]:
p1.add_incident(h_folder_1)

The provided incident name (hercules101) has been used before.
The incident name should be a unique name. Current incidents: ['hercules101', 'hercules102']

In tsprocess, the users are not allowed to directly process the records. However, they can define processing label. Each processing label has a unique name, a combination of the record and processing label generates a unique hash value which help us to track the processed records.

In [9]:
# let's take a look at valid processing labels
p1.valid_processing_labels()

0: lowpass_filter - args: {'fc': 'corner freq (Hz)', 'N': 'order (default:4)'}
1: highpass_filter - args: {'fc': 'corner freq (Hz)', 'N': 'order (default:4)'}
2: bandpass_filter - args: {'fcs': 'corner freqs (Hz)', 'N': 'order (default:4)'}
3: scale - args: {'factor': 'scaling factor'}
4: taper - args: {'flag': 'front, end, all', 'm': 'number of samples for tapering'}
5: cut - args: {'flag': 'front, end', 'm': 'number of samples for tapering', 't_diff': 'time to cut'}
6: zero_pad - args: {'flag': 'front, end', 'm': 'number of samples for tapering', 't_diff': 'time to add'}
7: rotate - args: {'angle': 'rotation angle'}
8: set_unit - args: {'unit': 'm, cm'}
9: set_vertical_or - args: {'ver_or': 'up, down'}
10: align_record - args: {'hc_or1': 'angle', 'hc_or2': 'angle', 'ver_or': 'up, down'}

Up to this point that I am preparing this report, there are 10 different processing labels and their types are self-explanatory. See API Reference page for more details.

In [10]:
# let's add some processing labels
# I want to add a lowpass filter with corner frequency at 2 Hz. Let's call it lpf2.
p1.add_processing_label('lpf2','lowpass_filter', {'fc':2, 'N':4})

Now, lpf2 is a reserved name and we cannot use it as a processing label name for other processing labels. For example, another attemp to reuse the name will be denied.

In [11]:
p1.add_processing_label('lpf2','highpass_filter', {'fc':0.1, 'N':4})

label name: 'lpf2' has been already used. Try another name.

tsprocess merge the records from different incidents based on their geographical locations. At this point any record that are closer than 10 m will be estimated at the same location. In the future versions, we will make it adjustable by the user. Choosing station for any process are based on stations filter. There are several approach to choose stations. See API Reference page for more details.

In [12]:
# let's take a look at valid stations filter (Not to be confused with timeseries filtering)
p1.valid_station_filter_type()

0: epi_dist_lt - args: {'distance': 'in km'}
1: epi_dist_gt - args: {'distance': 'in km'}
2: epi_dist_lte - args: {'distance': 'in km'}
3: epi_dist_gte - args: {'distance': 'in km'}
4: azimuth_bt - args: {'azmth': '[az1, az2]'}
5: include_stlist_by_incident - args: {'incident_name': 'name of incident', 'stations': 'list of stations'}
6: exclude_stlist_by_incident - args: {'incident_name': 'name of incident', 'stations': 'list of stations'}
In [13]:
# let's add some station filter. One can use several station filters.
# The final result will be the intersection of all used filters.
# For example, I want to study stations in less than 10 km from the source of earthquake.
# let's call it 'lesst10km'
p1.add_station_filter('lesst10km','epi_dist_lt',{'distance':10})

One can define station filter and processing label at any time of the processing (of course before using them). In this section, I define a couple of processing labels and station filters for future use.

In [14]:
p1.add_processing_label('stcm','set_unit',{'unit':'cm'}) # to convert all data to cm
p1.add_processing_label('hpf0.1','highpass_filter',{'fc':0.1, 'N':4}) # highpass filter with corner freq at 0.1 Hz.
p1.add_processing_label('bpf1_2','bandpass_filter',{'fcs':[1,2], 'N':4}) #bandpass filters with corner freqs 1,2 Hz.
p1.add_processing_label('taper_it','taper',{'flag':'all', 'm':100}) # tapers both side of the signal, 100 data samples each.

p1.add_station_filter('moret10km','epi_dist_gt',{'distance':10}) # stations at epicentral distance more than 10 km
p1.add_station_filter('lesst20km','epi_dist_lt',{'distance':20}) # stations at epicentral distance less than 20 km
p1.add_station_filter('azb0_90','azimuth_bt',{'azmth':[0,90]}) # stations at azimuth between 0 and 90 degrees.
p1.add_station_filter("select_stations", "include_stlist_by_incident",
                      {"incident_name":"hercules101",
                       "stations":['station.9','station.29','station.34']}) # select station location based on their names at one of incidents.
p1.add_station_filter("select_station", "include_stlist_by_incident",
                      {"incident_name":"hercules101",
                       "stations":['station.20']}) # select station location based on their names at one of incidents.

Now, we have enough materials to process data. At any time, that we need to add any processing labels or stations filters we can add them. let's take a look at those stations that are located between 10 and 20 km from the source of earthquake.

In [15]:
p1.show_stations_on_map2(['hercules101'],[[]],['moret10km','lesst20km'],{'llrtlatlon':[33.0,-119.5,35.5,-116.5]})

If we do not use any station filter, it will represent all available stations.

In [16]:
p1.show_stations_on_map2(['hercules101'],[[]],[],{'llrtlatlon':[33.0,-119.5,35.5,-116.5]})

We can take a look at the location of select stations.

In [17]:
p1.show_stations_on_map2(['hercules101'],[[]],['select_stations'],{'llrtlatlon':[33.0,-119.5,35.5,-116.5]})

We can take a look at the stations at azimuth between 0 and 90 degrees.

In [18]:
p1.show_stations_on_map2(['hercules101'],[[]],['azb0_90'],{'llrtlatlon':[33.0,-119.5,35.5,-116.5]})

Any processing with these station filters will be applied to those stations showed in the map.

In [19]:
# let's take a look at available processing labels.
p1.list_of_processing_labels()

lpf2 --> ['lowpass_filter', {'fc': 2, 'N': 4}]
hpf0.1 --> ['highpass_filter', {'fc': 0.1, 'N': 4}]
bpf1_2 --> ['bandpass_filter', {'fcs': [1, 2], 'N': 4}]
taper_it --> ['taper', {'flag': 'all', 'm': 100}]
stcm --> ['set_unit', {'unit': 'cm'}]
In [20]:
# let's take a look at velocity timeseries at select station (station.20) for 'hercules101'.
p1.plot_velocity_records(['hercules101'],[[]],['select_station'],{})
In [21]:
# let's have both of them ('hercules101' and 'hercules102')
p1.plot_velocity_records(['hercules101','hercules102'],[[],[]],['select_station'],{})
In [22]:
# This two simulation are done for one earthquake, but they seem very different. 
# Let's see what is the difference between these simulations.
import pandas as pd
table = p1.compare_incidents(['hercules101','hercules102'], only_differences=False)
with pd.option_context('display.max_rows', None, 'display.max_columns', None):  
    display(table)
parameters hercules101 hercules102
0 nonlinear_shear_velocity_cut 0 0
1 material_properties_type alphakay alphakay
2 do_damping_statistics 0 0
3 inputfiles_parameters inputfiles/parameters.in inputfiles/parameters.in
4 output_planes_print_rate 10 10
5 material_model linear linear
6 use_checkpoint 0 0
7 output_stations_directory outputfiles/stations outputfiles/stations
8 source_directory_output outputfiles/srctmp outputfiles/srctmp
9 print_station_velocities yes yes
10 incident_unit m m
11 simulation_delta_time_sec 0.001 0.001
12 region_length_east_m 28000.0 28000.0
13 region_origin_longitude_deg -118.20819 -118.20819
14 source_directory inputfiles/sourcefiles_extended inputfiles/sourcefiles_extended
15 material_properties_count 0 0
16 stiffness_calculation_method effective effective
17 simulation_end_time_sec 60.001 60.001
18 region_depth_deep_m 14000.0 14000.0
19 simulation_velocity_profile_freq_hz 0.0 0.0
20 mesh_coordinates_for_matlab no no
21 fixedbase_input_dir results/3rd-round/00bldgs/stations results/3rd-round/00bldgs/stations
22 drm_directory outputfiles/drm outputfiles/drm
23 implement_drm no no
24 print_station_accelerations yes yes
25 output_stats_file output-stats.txt output-stats.txt
26 geostatic_loading_time_sec 0 0
27 number_output_planes 1 1
28 output_planes_directory outputfiles/planes outputfiles/planes
29 softening_factor 4 4
30 checkpointing_rate 0 0
31 use_infinite_qk yes yes
32 number_of_buildings 0 0
33 fixedbase_input_startindex 000 000
34 output_parallel 0 0
35 source_hypocenter 33.933, -117.916, 5.1 33.933, -117.916, 5.1
36 checkpoint_path outputfiles/checkpoints outputfiles/checkpoints
37 hr_comp_orientation_2 129.9 129.9
38 output_velocity_file vel.h4d vel.h4d
39 which_drm_part part1 part1
40 hr_comp_orientation_1 39.9 39.9
41 output_displacement_file disp.h4d disp.h4d
42 drm_print_rate 0 0
43 simulation_displacement_out 0 0
44 fixedbase_input_dt 0.00 0.00
45 buildings_n_factor 0 0
46 output_stations_print_rate 10 10
47 output_velocity 0 0
48 include_nonlinear_analysis no no
49 the_threshold_Vp_over_Vs 3.0 3.0
50 number_output_stations 356 356
51 incident_source_hypocenter (33.933, -117.916, 5.1) (33.933, -117.916, 5.1)
52 drm_offset_x 0 0
53 region_depth_shallow_m 0.00000000 0.00000000
54 simulation_wave_max_freq_hz 4.0 4.0
55 simulation_output_rate 0 0
56 output_displacement 0 0
57 drm_offset_y 0 0
58 part1_delta_t 0.0 0.0
59 geostatic_cushion_time_sec 0 0
60 monitor_file outputfiles/monitor.txt outputfiles/monitor.txt
61 ver_comp_orientation down down
62 output_mesh 0 0
63 incident_name hercules101 hercules102
64 consider_fixed_base no no
65 simulation_velocity_out 0 0
66 planes_input_file inputfiles/parameters.in inputfiles/parameters.in
67 incident_folder hercules_1 hercules_2
68 mesh_coordinates_directory_for_matlab outputfiles/matlabmesh outputfiles/matlabmesh
69 simulation_shear_velocity_min 500 500
70 drm_edgesize 0 0
71 mesh_etree_output_file outputfiles/mesh.e outputfiles/mesh.e
72 the_threshold_damping 0.05 0.05
73 surface_shift_m 0 0
74 incident_type hercules hercules
75 fixedbase_input_sufix station station
76 region_azimuth_leftface_deg 0 0
77 4D_output_file outputfiles/disp.out outputfiles/disp.out
78 region_length_north_m 28000.0 28000.0
79 use_progressive_meshing 4 4
80 simulation_start_time_sec 0 0
81 include_buildings no no
82 region_origin_latitude_deg 33.85173 33.85173
83 simulation_node_per_wavelength 16 16
84 cvmdb_input_file /u/sciteam/khoshnev/scratch/eqsims/cvms/etree_... /u/sciteam/khoshnev/scratch/eqsims/cvms/etree_...
85 print_matrix_k no no
86 min_octant_size_m 0 0
87 number_profile_layers 1551 1551
88 type_of_damping none bkt
In [23]:
# this is too much and most of them are the same. 
# Just show the differences
p1.compare_incidents(['hercules101','hercules102'], only_differences=True)
Out[23]:
parameters hercules101 hercules102
0 incident_name hercules101 hercules102
1 incident_folder hercules_1 hercules_2
2 type_of_damping none bkt
In [24]:
# Ok, Now we can see the reason for very high amplitudes. Damping was not active at the first simulation.
# Let's take a look at the results, with applying lowpass filter (f=2)
p1.plot_velocity_records(['hercules101','hercules102'],[['lpf2'],['lpf2']],['select_station'],{})
In [25]:
# Now it makes sense, however, there is an artifact.
# let's first taper the signals then fitler them. 
p1.plot_velocity_records(['hercules101','hercules102'],
                         [['taper_it','lpf2'],['taper_it','lpf2']],
                         ['select_station'],
                         {})
In [26]:
# Can we take a closer look at the frequnecy content may be 0 - 4 ?
p1.plot_velocity_records(['hercules101','hercules102'],
                         [['taper_it','lpf2'],['taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_freq':[0,4]})
In [27]:
# How about time, can we take a look at between 5 and 12 seconds?
p1.plot_velocity_records(['hercules101','hercules102'],
                         [['taper_it','lpf2'],['taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_freq':[0,4], 'zoom_in_time':[5,12]})
In [28]:
# How is the acceleration timeseries look like?
p1.plot_acceleration_records(['hercules101','hercules102'],
                         [['taper_it','lpf2'],['taper_it','lpf2']],
                         ['select_station'],
                         {})
In [29]:
# I wonder to know how much filter order can affect the results?
# Especifically, is there any difference between N=4, N=8, and N=16? Let's say for lowpass filter at 2 Hz.
p1.add_processing_label('lpf2_n8','lowpass_filter',{'fc': 2, 'N': 8})
p1.add_processing_label('lpf2_n16','lowpass_filter',{'fc': 2, 'N': 16})
In [30]:
p1.plot_velocity_records(['hercules101','hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2_n8'],['taper_it','lpf2_n16']],
                         ['select_station'],
                         {'zoom_in_freq':[0,4]})

As we expected, higher filter order gives sharp cut at the corner frequency.

In [31]:
# And by the way, how is this station components orientation?
p1.plot_records_orientation(['hercules101','hercules102'],
                         [[],[]],
                         ['select_station'],
                         {})
In [32]:
# Intersting! can we switch the vertical component to up?
p1.add_processing_label('v_up','set_vertical_or',{'ver_or':'up'})
In [33]:
p1.plot_records_orientation(['hercules101','hercules102'],
                         [[],['v_up']],
                         ['select_station'],
                         {})
In [34]:
# I do not beleive you :). Can I see the timeseries?
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['v_up','taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_freq':[0,4], 'zoom_in_time':[0,20]})
In [35]:
# That is right. 
# Can I align the record to 0 and 90 degrees?
p1.add_processing_label('align_0_90','align_record',{'hc_or1':0, 'hc_or2':90, 'ver_or':'down'})
In [36]:
p1.plot_records_orientation(['hercules101','hercules101'],
                         [[],['align_0_90']],
                         ['select_station'],
                         {})
In [37]:
# Great! I assume, 0,90, down and 180, 270, up should be in the upozit direction, right?
# let's try:
p1.add_processing_label('align_180_270','align_record',{'hc_or1':180, 'hc_or2':270, 'ver_or':'up'})
In [38]:
p1.plot_records_orientation(['hercules101','hercules101'],
                         [['align_0_90'],['align_180_270']],
                         ['select_station'],
                         {})
In [39]:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['align_0_90','taper_it','lpf2'],['align_180_270','taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_time':[1,20],'zoom_in_freq':[0,4]})
In [40]:
# OK. Cool. Can we take a look at the same comparison, but on the select stations.
p1.plot_velocity_records(['hercules102','hercules102'],
                         [['align_0_90','taper_it','lpf2'],['align_180_270','taper_it','lpf2']],
                         ['select_stations'],
                         {'zoom_in_time':[1,20], 'zoom_in_freq':[0,4]})

Here are some other processing that is possible to conduct on time series.

In [41]:
# Scaling
p1.add_processing_label('sc1.2','scale',{'factor':1.2})
In [42]:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2','sc1.2'],['taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_time':[1,20],'zoom_in_freq':[0,4]})
In [43]:
# Zero-padding at the fron for one second.
p1.add_processing_label('zp1.0f','zero_pad',{'flag':'front', 'm':100, 't_diff':1})
In [44]:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2','zp1.0f'],['taper_it','lpf2']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [45]:
# rotate time series (45 degrees) 
p1.add_processing_label('rt45','rotate',{'angle':45})
In [46]:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2','rt45']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [47]:
# Let's take a look at orientation 
p1.plot_records_orientation(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2','rt45']],
                         ['select_station'],
                         {})
In [48]:
# Let's define some other rotation angles.

p1.add_processing_label('rt10','rotate',{'angle':10})
p1.add_processing_label('rt90','rotate',{'angle':90})
p1.add_processing_label('rt180','rotate',{'angle':180})
p1.add_processing_label('rt360','rotate',{'angle':360})
In [49]:
# So, if I rotate the record 8 times by 45 degrees I should get the same results.
# let's see:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2','rt45','rt45','rt45','rt45','rt45','rt45','rt45','rt45']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [50]:
# So, if I rotate the record 4 times by 90 degreess we should get the same result.
# let's see:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2','rt90','rt90','rt90','rt90']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [51]:
# So, results of rotating the signla by 90 degrees, and twice by 45 degrees and 9 times by 10 degrees should be the same.
# let's see:
p1.plot_velocity_records(['hercules101','hercules101','hercules101'],
                         [['taper_it','lpf2','rt90'],
                          ['taper_it','lpf2','rt45','rt45'],
                          ['taper_it','lpf2','rt10','rt10','rt10','rt10','rt10','rt10','rt10','rt10','rt10']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [52]:
p1.plot_records_orientation(['hercules101','hercules101','hercules101','hercules101','hercules101','hercules101'],
                         [['taper_it','lpf2','rt10'],
                          ['taper_it','lpf2','rt10','rt10'],
                          ['taper_it','lpf2','rt10','rt10','rt10','rt10'],
                          ['taper_it','lpf2','rt10','rt10','rt10','rt10','rt10'],
                          ['taper_it','lpf2','rt10','rt10','rt10','rt10','rt10','rt10'],
                          ['taper_it','lpf2','rt10','rt10','rt10','rt10','rt10','rt10','rt10'],
                         ],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [53]:
p1.plot_velocity_records(['hercules101','hercules101'],
                         [['taper_it','lpf2'],['taper_it','lpf2','rt360']],
                         ['select_station'],
                         {'zoom_in_time':[0,30],'zoom_in_freq':[0,4]})
In [54]:
# We conducted many analyses. Are they stored in the database?
# Yes. We can take a look at a summary of the database:
p1.database_summary()

Database name: project_1_db.sqlite
Database size: 400.285696 MB.
Number of items: 415
List of database trackers: ['project_1_dbtracker']
In [55]:
# I think this is enough for example 1. We close the database and close this page.
# Closing database makes sure that there is no pending data in the tracking buffer.
p1.close_database()

Database (project_1_db.sqlite) is closed.