How to automatically enquire the availability of seismic data using Obspy

Utpal Kumar   2 minute read   
   visitor badge  

In this post, we will see how can we retrive the available seismic waveforms information for a given network, station, channel and client in a given period of time.

This post will see how we can retrieve the available seismic waveforms information for a given network, station, channel, and client in a given period.

Retrieve all the clients’ list from Obspy

First, we need to retrieve the available clients from Obspy. This we can do using the obspy.clients.fdsn.header module.

from obspy.clients.fdsn.header import URL_MAPPINGS
for key in sorted(URL_MAPPINGS.keys()):
    print("{0:<11} {1}".format(key,  URL_MAPPINGS[key]))
all_clients = list(URL_MAPPINGS.keys())
all_clients.remove('IRIS')
all_clients = ['IRIS'] + all_clients

In the above code, we first retrieve the available clients from the Obspy fdsn service. This list is alphabetic, but we want to give preference to the “IRIS” client. So, I made the “IRIS” as the first element in the list.

The available list of the clients at the time of this post is:

IRIS        http://service.iris.edu
ISC         http://isc-mirror.iris.washington.edu
KNMI        http://rdsa.knmi.nl
KOERI       http://eida.koeri.boun.edu.tr
LMU         http://erde.geophysik.uni-muenchen.de
NCEDC       http://service.ncedc.org
NIEP        http://eida-sc3.infp.ro
NOA         http://eida.gein.noa.gr
ODC         http://www.orfeus-eu.org
ORFEUS      http://www.orfeus-eu.org
RASPISHAKE  http://fdsnws.raspberryshakedata.com
RESIF       http://ws.resif.fr
SCEDC       http://service.scedc.caltech.edu
TEXNET      http://rtserve.beg.utexas.edu
UIB-NORSAR  http://eida.geo.uib.no
USGS        http://earthquake.usgs.gov
USP         http://sismo.iag.usp.br

Retrieve the station information

# Define parameters
starttime = UTCDateTime("2020-05-15T00:00:00Z")
endtime = starttime+7*24*3600
net = "NC"
stn = "*"
channel = "*Z"
count = 0
success = False
for cl in all_clients:
    try:
        print(f"--> Trying for client: {cl}")
        client = Client(cl)

        inventory = client.get_stations(network=net, station=stn, channel=channel,
                                        level="response", starttime=starttime, endtime=endtime)

        print(inventory)

        inventory.write('station_info.txt', 'STATIONTXT', level='station')
        success = True
        if success:
            try:
                plot_stations()
            except:
                sys.exit()
            break
    except KeyboardInterrupt:
        sys.exit()
    except:
        print(cl, sys.exc_info())
    count += 1

The above code will try to retrieve the station information up to the “response” level from the list of clients iteratively. If it successfully obtains the information, it will break the loop and plot the stations on a map.

I enquired for all stations in the “NC” network and vertical component. The time range of the waveform is seven days from “2020-05-15”. All these parameters can be modified depending on the needs. The obtained station inventory will be written in the file “station_info.txt”.

Plot the stations

Now, we read the “station_info.txt” as a pandas dataframe and plot the stations using the “pygmt”. If multiple networks are enquired, then they will be plotted in a different color (colors are randomly chosen).

import numpy as np
import pygmt
import pandas as pd
np.random.seed(45)  # to get the same color at each run


def plot_stations():
    df = pd.read_csv('station_info.txt', delimiter='|')
    print(df.head())

    # get the list of networks
    networks = list(set(df['#Network'].tolist()))

    dfs = []
    for net in networks:
        df1 = df[df['#Network'] == net]
        dfs.append(df1)

    colorsList = []
    for i in range(len(networks)):
        colorsList.append('#%06X' % np.random.randint(0, 0xFFFFFF))

    minlon, maxlon = df['Longitude'].min()-1, df['Longitude'].max()+1
    minlat, maxlat = df['Latitude'].min()-1, df['Latitude'].max()+1

    # define etopo data file

    topo_data = "@earth_relief_30s"

    # Visualization
    fig = pygmt.Figure()
    # make color pallets
    pygmt.makecpt(
        cmap='etopo1',
        series='-8000/8000/1000',
        continuous=True
    )

    # plot high res topography
    fig.grdimage(
        grid=topo_data,
        region=[minlon, maxlon, minlat, maxlat],
        projection='M4i',
        shading=True,
        frame=True
    )

    # plot coastlines
    fig.coast(
        region=[minlon, maxlon, minlat, maxlat],
        projection='M4i',
        shorelines=True,
        frame=True
    )

    for idx, dff in enumerate(dfs):
        fig.plot(
            x=dff["Longitude"].values,
            y=dff["Latitude"].values,
            style="i10p",
            color=colorsList[idx],
            pen="black",
            label=networks[idx]
        )

    fig.legend(position="JTR+jTR+o0.2c", box=True)

    fig.savefig('station_map.png', crop=True, dpi=300)


if __name__ == '__main__':
    plot_stations()

Next, if you want to download the data, you can check my post for downloading seismic waveforms.

Plot of the stations with data availability
Plot of the stations with data availability

Complete codes

The complete code for this work can be downloaded from my github repository.

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.


Leave a comment