Getting started with ObsPy

#1

I’m in the process of learning to use Python and ObsPy on a Win10 machine. I have almost no experience with Python, so I have worked through several tutorials to get acquainted with the basics.

I went to Anaconda to get ObsPy. Through Anaconda I installed ObsPy (I think). However, I am just not able to access the library – operator error I’m sure.

At this point, I really don’t know what I’m doing! I’d like to reach out to the Shakeosphere and get some startup support using ObsPy on a Windows machine.

Once off and running, I plan to process the raw seismic files for plotting, filtering, etc. Basically just tinkering, but I need to get out of the box first!

Thanks in advance,
Jeff

#2

Hi Jeff, I’m psyched that you’re interested in getting your feet wet in computational seismology!

Did you follow the instructions on the obspy website for installing via Anaconda? https://github.com/obspy/obspy/wiki/Installation-via-Anaconda

What happens when you open an Anaconda Prompt and paste and run the following commands?

conda activate obspy
python -c 'import obspy; print(obspy.__version__)'
#3

The activate command changes the prompt from (base) to (obspy), so that seems to work. The request to view version # returns with ModuleNotFoundError:

(base) PS C:\Users\Jeff> conda activate obspy
(obspy) PS C:\Users\Jeff> python -c ‘import obspy; print(obspy.version)’
Traceback (most recent call last):
File “”, line 1, in
ModuleNotFoundError: No module named ‘obspy’
(obspy) PS C:\Users\Jeff>

#4

Here are a few screen shots for reference. The first is the Anaconda Powershell where I ran the commands. The second shows the directory where obspy was installed by Anaconda, and the last shows Anaconda Navigator with reference to obspy:

#5

Ok it looks like you have the environment but obspy is not installed in it yet. Try this:

conda activate obspy
conda config --add channels conda-forge
conda install obspy
python -c 'import obspy; print(obspy.__version__)'
#6

Success! Version 1.1.1 returned.
OK, what did I just do?
What do I have to do to get started each time?

#7

Here’s my first obspy script.
It reads, filters and plots a single miniSEED file from the RS1D.
Apologies for the display here, all my text gets interpreted.
At the bottom of this message is an example plot of 24 hrs recording with a 5 Hz lo pass!!!

Any suggestions/upgrades welcome!
(edit: some edits suggested by @iannesbitt below)


# Plot and filter miniSEED

# Activate obspy from root:

     # conda activate obspy

# Note: there's actually no need to do "import obspy" here because
#    you only use the read() function within obspy. So you can save
#    a little bit of computing time and memory by not importing the
#    whole thing!
from obspy import read

# Read file and print stats

st = read("D:\JLL Data Files\PDF_files\EQ Software\Shake R2DEF miniSEED files\AM.R2DEF.00.EHZ.D.2019.214")
tr = st[0]
msg =  "%s %s %s %f %f" % (tr.stats.station, str(tr.stats.starttime), str(tr.stats.endtime), tr.data.mean(), tr.data.std())
print("\n\n",msg)

# Filter and Plot

st.filter("lowpass", freq=5.0)
st.plot()

2 Likes
#8

Great! Here’s an explanation:

# this activates the obspy environment you created earlier.
#    now that you have obspy installed, this is the only thing you
#    need to do before running your script.
conda activate obspy

# this is just a one-time command that adds the conda-forge channel
#    to the places that conda searches for software in.
#    conda-forge is where obspy lives.
conda config --add channels conda-forge

# this (obviously) installs obspy into the obspy environment
conda install obspy

# this runs a command (-c) quickly in Python without the need for you
#    to fully dive into a Python shell. so it opens a python shell, runs
#    anything inside the single quotes, then closes the shell again.
#    the two commands (separated by a semicolon) actually import obspy
#    and print the version, but only inside that briefly existing shell.
python -c 'import obspy; print(obspy.__version__)'

Given that, here’s how you start up and run your script:

  1. Open an Anaconda Powershell
  2. Activate the (obspy) environment
    conda activate obspy
    
  3. Run your script:
    python myscript.py
    

Simple, right?

#9

I also made some formatting and code suggestion edits to your post above. Hope these are helpful and clarifying.

#10

Moving along!

Here I access waveform and instrument response from Raspberryshakedata and apply (remove) the instrument response.

I am confused about the results. The RS1D geophone should be in velocity units, but when I ask for VEL, the seismogram is riding on a long period signal. When I ask for ACC, I get a seismogram that is higher frequency and noisier than the original.

Below are the plots: Black curve (1) = input, Green curve (2) = ACC, Blue curve (3) = VEL, and Red curve (4)= DISP.

Here is the script:

# remove instrument response

# Load stuff
from obspy import read_inventory, read
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream

# set the station name and download the response information
stn = 'R2DEF'
inv = read_inventory('https://fdsnws.raspberryshakedata.com/fdsnws/station/1/query?network=AM&station=%s&level=resp&format=xml' % (stn))
print (inv)

# set data start/end times
start = UTCDateTime(2019, 8, 1, 14, 25, 30) # (YYYY, m, d, H, M, S)
end = UTCDateTime(2019, 8, 1, 14, 27, 0) # (YYYY, m, d, H, M, S)

# set the FDSN server location and channel names
rs = Client(base_url='https://fdsnws.raspberryshakedata.com/')
channels = ['EHZ'] # ENx = accelerometer channels; EHx or SHZ = geophone channels

# get waveforms, put them all into one Stream, and attach the response
stream = Stream()
for ch in channels:
    trace = rs.get_waveforms('AM', stn, '00', ch, start, end)
    stream += trace
stream.plot(color='black', linewidth='0.5',handle='true')

# attach the response info
stream.attach_response(inv)

# remove the response and plot ACC
resp_removed_ACC = stream.remove_response(output='ACC') # convert to acceleration in M/S
resp_removed_ACC.plot(color='darkgreen', linewidth='0.5', handle='true')

# remove the response and plot VEL
resp_removed_VEL = stream.remove_response(output='VEL') # convert to acceleration in M/S
resp_removed_VEL.plot(color='darkblue', linewidth='0.5', handle='true')

# remove the response and plot DISP
resp_removed_DISP = stream.remove_response(output='DISP') # convert to acceleration in M/S
resp_removed_DISP.plot(color='darkred', linewidth='0.5')
1 Like
#11

Hi @OBSeismo, this is quickly reaching the limits of my expertise area. Your ACC plot looks okay to me but can you try plotting DISP and VEL after highpass filtering above 0.5 Hz? It looks like the long-period stuff is an artifact of some kind (although comparing it to a nearby broadband might be informative).