Section Plot error

Usage Questions
1

Dear Obspy forum,

I have been trying to plot a “section plot” with data that comes from Nordic format (Seisan), at least to plot the vertical (Z) and one horizontal (N or E) channels, the Nordic file has 15 channels from 5 seismic stations, so to try to do this I have made this code:


import os
from obspy.core import read
from obspy.core import AttribDict
from obspy import UTCDateTime
from obspy import Trace


os.chdir('/home/asus/Documents/Pdoc_res/PhN_picks/ML1.0')

trace = Trace()

# Station SODC... coordinates

# Set the latitude and longitude for the station number 2, closest to epicenter
trace.stats.coordinates = {'latitude': -17.2057, 'longitude': -65.8726}

# Earthquake coordinates
latitude = -17.219
longitude = -65.858
ev_coord = {'latitude': latitude, 'longitude': longitude}

# Read the waveform data from a Nordic format file
stream = read("2022-05-18-2002-42S.NSN___010")

# Add P and S phase arrival times from the station #2 (closest to epicenter)
p_arrival = UTCDateTime(2022, 5, 18, 20, 2, 50.2)  #  P arrival time
s_arrival = UTCDateTime(2022, 5, 18, 20, 2, 52.4)  # S arrival time

for trace in stream:
    trace.stats.sac = AttribDict()  # Create a SAC header
    trace.stats.sac.t1 = p_arrival - trace.stats.starttime
    trace.stats.sac.t2 = s_arrival - trace.stats.starttime

# Max distance from epicenter to stations
distance_km = 15
distance_m = distance_km * 1000
trace.stats.distance = distance_m

# Try to plot the waveform data with P and S phase markers
stream.plot(type='section', show=True, draw_neutral=False, 
            time_down=True, offset_min=-10, offset_max=10, dist_degree=True,
            ev_coords=(ev_coord)

However, I obtain this error:

ValueError: Define latitude/longitude in trace.stats.coordinates and ev_coord. See documentation.

I looked the docs but I don’t understand where is my error, would you mind please tell me what did I do wrong?
Thanks in advance.
2022-05-18-2002-42S.NSN___010 (122.2 KB)

2

It looks like you have not defined the trace.stats.coordinate.latitude and trace.stats.coordinate.longitude for your given stream. You have also passed ev_coords as a tuple containing a dictionary. It isn’t very clear to me what you are trying to do with the sac headers or setting the distance, so apologies if I have misunderstood.

In general for this kind of a plot you would want to have some way of associating a trace to a given location (in the below example I use a dictionary to make this explicit, but it would be simple with an Inventory if you have one), and the location of your source.

Example (not tested and will need to be modified for your locations):

from obspy import read

# To fill in for your stations, a dictionary of tuples of (latitude, longitude) keyed 
# by seed ID without channel. Alternatively you could use an Inventory and use 
# the .select method to get the correct locations
station_locations = {
    "net.sta.loc": (-45.0, 170.0),
}

st = read("your_data")

for tr in st:
    # If you get a KeyError here you probably don't have the station in your station_locations dict
    loc = station_locations[f"{tr.stats.network}.{tr.stats.station}.{tr.stats.location}"]
    tr.stats.coordinates.latitude = loc[0]
    tr.stats.coordinates.longitude = loc[1]

# To complete with your event latitude and longitude
ev_coord = (ev_lat, ev_long)

st.plot(type="section", ev_coord=ev_coord, dist_degree=True)

If you wanted to plot your record section in meters from the source you would need to calculate the distance from your source to your station for each station and set tr.stats.distance for each trace to the correct calculated value. You can use obspy.geodetics.gps2dist_azimuth to calculate great circle distances.

1 Like
3

@calum-chamberlain , thank you. I thought to add P and S time to waveform header’s to sort the waveforms, for that reason I tried to use SAC format.
I will try your code and let you know.
Stay safe and best regards,
Tonino