GPS (GNSS) data handling

Loading a GNSS campaign from CSV files

In this example we load GPS station locations and displacement data from the 1994 Northridge Earthquake into a GNSSCampaign data model. The model is then stored in a YAML file and loaded again using pyrocko.guts.

Download gnss_campaign.py

from pyrocko import guts
from pyrocko.model import gnss
from pyrocko.example import get_example_data

# Data from Hudnut et al. 1996

# K. W. Hudnut, Z. Shen, M. Murray, S. McClusky, R. King, T. Herring,
# B. Hager, Y. Feng, P. Fang, A. Donnellan, Y. Bock;
# Co-seismic displacements of the 1994 Northridge, California, earthquake.
# Bulletin of the Seismological Society of America; 86 (1B): S19-S36. doi:

# https://pasadena.wr.usgs.gov/office/hudnut/hudnut/nr_bssa.html

mm = 1e3

fn_stations = 'GPS_Stations-Northridge-1994_Hudnut.csv'
fn_displacements = 'GPS_Data-Northridge-1994_Hudnut.csv'

get_example_data(fn_stations)
get_example_data(fn_displacements)

campaign = gnss.GNSSCampaign(name='Northridge 1994')

# Load the stations
with open(fn_stations, 'r') as f:
    for line in f:
        if line.startswith('#'):
            continue
        row = line.split(',')
        sta = gnss.GNSSStation(
            code=row[0].strip(),
            lat=float(row[1]),
            lon=float(row[2]),
            elevation=float(row[3]))
        campaign.add_station(sta)

# Load the displacements
with open(fn_displacements, 'r') as f:
    for line in f:
        if line.startswith('#'):
            continue
        row = line.split(',')

        station_id = row[0].strip()
        station = campaign.get_station(station_id)

        station.east = gnss.GNSSComponent(
            shift=float(row[1]) / mm,
            sigma=float(row[2]) / mm)
        station.north = gnss.GNSSComponent(
            shift=float(row[7]) / mm,
            sigma=float(row[8]) / mm)
        station.up = gnss.GNSSComponent(
            shift=float(row[14]) / mm,
            sigma=float(row[15]) / mm)

print('Campaign %s has %d stations' % (campaign.name, campaign.nstations))
campaign.dump(filename='GPS_Northridge-1994_Hudnut.yml')

# Load the campaign back in
campaign_loaded = guts.load(filename='GPS_Northridge-1994_Hudnut.yml')

Loading and mapping of GNSS campaign data from UNR

The Nevada Geodetic Laboratory (UNR) releases co-seismic GNSS surface displacements for significant earthquakes. This script shows the import of such co-seismic displacement tables for the 2019 Ridgecrest earthquake and mapping through automap.

Download gnss_unr_campaign.py

#!/usr/bin/env python3
import numpy as num
import requests

from pyrocko.model import gnss, location
from pyrocko.plot.automap import Map

km = 1e3
DIST_MAX = 80. * km
RIDGECREST_EQ = location.Location(lat=35.766, lon=-117.605)
GNSS_URL = 'http://geodesy.unr.edu/news_items/20190707/ci38457511_forweb.txt'

fname = 'ci38457511_forweb.txt'
# Download the data
with open(fname, 'wb') as f:
    print('Downloading data from UNR...')
    r = requests.get(GNSS_URL)
    f.write(r.content)

campaign = gnss.GNSSCampaign(name='Ridgecrest coseismic (UNR)')
print('Loading Ridgecrest GNSS data from %s (max_dist = %.1f km)' %
      (fname, DIST_MAX / km))

# load station names and data
with open(fname, 'r') as f:
    for header in range(2):
        next(f)
    names = [l.split(' ')[0] for l in f]
gnss_data = num.loadtxt(fname, skiprows=2, usecols=(1, 2, 3, 4, 5, 6, 7, 8))

for ista, sta_data in enumerate(gnss_data):
    name = names[ista]
    lon, lat, de, dn, du, sde, sdn, sdu = map(float, sta_data)

    station = gnss.GNSSStation(
        code=name,
        lat=lat,
        lon=lon,
        elevation=0.)

    if station.distance_to(RIDGECREST_EQ) > DIST_MAX:
        continue

    station.east = gnss.GNSSComponent(
        shift=de,
        sigma=sde)
    station.north = gnss.GNSSComponent(
        shift=dn,
        sigma=sdn)
    station.up = gnss.GNSSComponent(
        shift=du,
        sigma=sdu)

    campaign.add_station(station)

print('Selected %d stations, saving to ridgecrest_coseis_gnss.yml'
      % (campaign.nstations))
campaign.dump(filename='2019-ridgecrest_gnss_unr.yml')

print('Creating map...')
m = Map(
    lat=RIDGECREST_EQ.lat,
    lon=RIDGECREST_EQ.lon,
    radius=DIST_MAX * 1.5,
    width=20.,
    height=14.)
m.add_gnss_campaign(campaign)
m.save('2019-ridgecrest_GNSS-UNR.pdf')