Comparing results from response removal with results from SAC software

I am trying to download some seismic data in SAC format. I have used the following python code for the download.

# Loop over the events
for event in events:
      # Get event information
      event_time =[0].time
      event_latitude =[0].latitude
      event_longitude =[0].longitude
      event_magnitude = event.magnitudes[0].mag
      event_time_str = event_time.strftime('%Y%m%d%H%M%S')

    # Loop over available networks
    for network in event_client.get_stations(network=network_code, starttime=start_time, endtime=end_time, minlatitude=min_lat, maxlatitude=max_lat, minlongitude=min_lon, maxlongitude=max_lon):
        for station in network:
            station_code = station.code
            station_latitude = station.latitude
            station_longitude = station.longitude

            # Calculate the back azimuth, azimuth, distance, delta, and gcarc
            distance, azimuth, back_azimuth = gps2dist_azimuth(
                event_latitude, event_longitude, station_latitude, station_longitude
            distance = distance/1000
            gcarc = locations2degrees(
                event_latitude, event_longitude, station_latitude, station_longitude

            if not (3330 <= distance <= 11132):
                print(f"Distance between {station_code} and event {event_time} not in required range, skipping...")

                # Download waveform data for the station
                st = waveform_client.get_waveforms(
                    starttime=event_time - 300,  # Adjust for the desired time window
                    endtime=event_time + 1500,    # Adjust for the desired time window
            except obspy.clients.fdsn.header.FDSNNoDataException:
                print(f"No data available on {event_time} for station {station_code}, skipping...")
                continue  # Skip this station and move on to the next one

            # Check if the Stream is empty or has no traces
            if not st or len(st) == 0:
                print(f"No valid data for station {station_code}, skipping...")
                continue  # Skip this station and move on to the next one

            # Attempt to remove the response
            except Exception as e:
            # Handle the exception (e.g., skip the trace)
              print(f"Skipping trace for station {station_code}: {e}")

            # Create a directory for each event
            event_directory = os.path.join(output_directory, event_time_str)
            os.makedirs(event_directory, exist_ok=True)

            print(f"Data available on {event_time} for station {station_code}, proceeding...")

            st.taper(max_percentage=0.05, type='hann')
            st.filter('bandpass', freqmin=0.005, freqmax=1)

            # Save the processed data as SAC files in the event-specific folder
            for trace in st:
                network =
                station = trace.stats.station
                channel =
                location = trace.stats.location
                sample_rate = trace.stats.sampling_rate

                sac_trace = SACTrace.from_obspy_trace(trace)
       = azimuth
                sac_trace.baz = back_azimuth
                sac_trace.dist = distance
       = 1/sample_rate
                sac_trace.gcarc = gcarc
                sac_trace.mag = event_magnitude  # Set event magnitude
                sac_trace.o = 300
                sac_trace.stla = station_latitude  # Set station latitude
                sac_trace.stlo = station_longitude  # Set station longitude
                sac_trace.evla = event_latitude  # Set event latitude
                sac_trace.evlo = event_longitude  # Set event longitude

                sac_file_path = os.path.join(
                sac_trace.write(sac_file_path, byteorder="little")

            print(f"Data on {event_time} for station {station_code} saved, proceeding...")

# Print a message indicating the folder where all the data is saved
print(f"All data saved in the folder: {output_directory}")

This is just the main part of the code. I understand that removing response in VEL results in m/s as the units. But the values are in the 1e-16 range. Is this correct? This data is then used in some program do_mft to extract amplitude vs time period data. I am self-taught in Python and ObsPy. Any and all help is appreciated. I have attached an example plot below.

Looks like you are calling st.remove_response() twice, so your instrument response gets removed twice which isn’t what you want. That function like many others works in place to save memory, which is also stated in docs.

I’ve updated the pinned “Welcome” post with instructions for Markdown, please check that out to make code blocks render properly. :slight_smile:

Thanks a lot. This has helped. Is this more in line with what the values should be?

Are there any more points to be taken care off? Also if there any good guides to working with data in seismology, it would be helpful.

Upon comparing with some other sac files that i had received earlier, I found that the amplitude ranges are in the 1e3 to 1e4 magnitude range and in cm/s units. So why are these values much lower? Both are of similar magnitudes.

Please read the docs, it states that for velocity output is in SI units, i.e. in m/s. This has also been asked before in similar posts on the topic of instrument removal.

No. My doubt is are these values right. If I convert these values to cm/s the amplitude values lie around 1e-3 to 1e-4 range, while the sac files that I received have amplitude in cm/s and the values are around 1e4. This is why I have my doubts.

Instrument correction is not a binary “done/not done” operation. There are quite a few choices to be made in the details that can drastically change your output, you should especially look at the parameters water_level and pre_filt in remove_response() and there’s the plot=True option to get an idea of the various steps in response removal. Also you can look at sacsim parameter, I haven’t used that in ages but I believe we used that at the very early times of obspy when comparing response removal results to SAC (it should only change results very minorly, just needed to get 1:1 the same results, some changes in tapering and demeaning). Then I always recommend having a buffer around your time window of interest that is cut off after processing, because you will/can get artifacts at start and end.
There is also a very very old test case comparing response removal with results from SAC, that was before the high level function remove_response was implemented, so it’s using underlying routines, but you can look at those details if you want.

Lastly, we get a lot of questions around response removal all the time, please also use the (advanced) search to look at previously answered questions/threads with similar questions asked, thanks!
E.g.: Search results for 'response #usage-questions after:2022-01-01' - ObsPy Forum

Edit: Oh and thanks for fixing the markdown in your original post with the syntax highlighted code block… appreciated! :heart: