View GitHub Repository rsudp

Schedule Inspection
Subscribe
Passed
Pull Request — master (#64)
by
unknown
06:49
created

build.amplitude_alert.PGVAlert.process()   A

↳ Parent: build.amplitude_alert

Complexity

Conditions 2

Size

Total Lines 20
Code Lines 11

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 11
dl 0
loc 20
rs 9.85
c 0
b 0
f 0
cc 2
nop 2
1 
import obspy
2 
import numpy as np
3 
from rsudp.client import Client
4 

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    5
                

                                    
                                                    
                
                
                
class PGVAlert:

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    6
                

                                    
                                                    
                
                
                
    def __init__(self, threshold):

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    7
                

                                    
                                                    
                
                
                
        self.threshold = threshold  # Richter magnitude threshold

            


            
                                                                                                            
                                                                
            
                                    
            
            

                
                    8
                

                                    
                                                    
                
                
                

            


            
                                                                        
                            
            
                                    
            
            

                
                    9
                

                                    
                                                    
                
                
                
    def process(self, data):

            


            
                                                                        
                            
            
                                    
            
            

                
                    10
                

                                    
                                                    
                
                
                
        # Create a Stream object from the incoming data

            


            
                                                                        
                            
            
                                    
            
            

                
                    11
                

                                    
                                                    
                
                
                
        st = obspy.Stream(obspy.Trace(data))

            


            
                                                                        
                            
            
                                    
            
            

                
                    12
                

                                    
                                                    
                
                
                
        # Remove mean and linear trends

            


            
                                                                        
                            
            
                                    
            
            

                
                    13
                

                                    
                                                    
                
                
                
        st.detrend("demean")

            


            
                                                                        
                            
            
                                    
            
            

                
                    14
                

                                    
                                                    
                
                
                
        st.detrend("linear")

            


            
                                                                        
                            
            
                                    
            
            

                
                    15
                

                                    
                                                    
                
                
                
        # Apply bandpass filter to isolate frequencies of interest

            


            
                                                                        
                            
            
                                    
            
            

                
                    16
                

                                    
                                                    
                
                
                
        st.filter("bandpass", freqmin=0.1, freqmax=20.0)

            


            
                                                                        
                            
            
                                    
            
            

                
                    17
                

                                    
                                                    
                
                
                
        # Integrate to convert displacement to velocity

            


            
                                                                        
                            
            
                                    
            
            

                
                    18
                

                                    
                                                    
                
                
                
        st.integrate()

            


            
                                                                        
                            
            
                                    
            
            

                
                    19
                

                                    
                                                    
                
                
                
        # Find Peak Ground Velocity (PGV)

            


            
                                                                        
                            
            
                                    
            
            

                
                    20
                

                                    
                                                    
                
                
                
        pgv = max(abs(st[0].data))

            


            
                                                                        
                            
            
                                    
            
            

                
                    21
                

                                    
                                                    
                
                
                
        # Convert PGV from m/s to cm/s

            


            
                                                                        
                            
            
                                    
            
            

                
                    22
                

                                    
                                                    
                
                
                
        pgv_cm_s = pgv * 100

            


            
                                                                        
                            
            
                                    
            
            

                
                    23
                

                                    
                                                    
                
                
                
        # Estimate Richter magnitude using the empirical relationship

            


            
                                                                        
                            
            
                                    
            
            

                
                    24
                

                                    
                                                    
                
                
                
        magnitude = np.log10(pgv_cm_s) + 2.4

            


            
                                                                        
                            
            
                                    
            
            

                
                    25
                

                                    
                                                    
                
                
                

            


            
                                                                        
                            
            
                                    
            
            

                
                    26
                

                                    
                                                    
                
                
                
        # Trigger alert if magnitude exceeds threshold

            


            
                                                                        
                            
            
                                    
            
            

                
                    27
                

                                    
                                                    
                
                
                
        if magnitude >= self.threshold:

            


            
                                                                        
                            
            
                                    
            
            

                
                    28
                

                                    
                                                    
                
                
                
            self.trigger_alert(magnitude)

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    29
                

                                    
                                                    
                
                
                

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    30
                

                                    
                                                    
                
                
                
    def trigger_alert(self, magnitude):

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    31
                

                                    
                                                    
                
                
                
        print(f"Earthquake detected! Estimated Richter Magnitude: {magnitude}")

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    32
                

                                    
                                                    
                
                
                

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    33
                

                                    
                                                    
                
                
                
# Initialize and start the client with the custom alert

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    34
                

                                    
                                                    
                
                
                
client = Client()

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    35
                

                                    
                                                    
                
                
                
alert = PGVAlert(threshold=4.0)

            


            
                                                                                                            
                            
            
                                    
            
            

                
                    36
                

                                    
                                                    
                
                
                
client.add_listener(alert.process)

            


            
                                                                                                            
                                                                
            
                                    
            
            

                
                    37
                

                                    
                                                    
                
                
                
client.run()

            


            
                                                        
            
                                    
            
            

                
                    38