|
# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
import logging
import math
import numpy as num
from collections import OrderedDict
import pyrocko.orthodrome as od
from pyrocko.guts import (Object, Float, String, List, StringChoice,
DateTimestamp)
from pyrocko.model import Location
guts_prefix = 'pf.gnss'
logger = logging.getLogger('pyrocko.model.gnss')
class GNSSComponent(Object):
''' Component of a GNSSStation
'''
unit = StringChoice.T(
choices=['mm', 'cm', 'm'],
help='Unit of displacement',
default='m')
shift = Float.T(
default=0.,
help='Component\'s shift in unit')
sigma = Float.T(
default=0.,
help='One sigma uncertainty of the measurement')
def __add__(self, other):
if not isinstance(other, self.__class__):
raise AttributeError('Other has to be of instance %s'
% self.__class__)
comp = self.__class__()
comp.shift = self.shift + other.shift
comp.sigma = math.sqrt(self.sigma**2 + other.sigma**2)
return comp
def __iadd__(self, other):
self.shift += other.shift
self.sigma = math.sqrt(self.sigma**2 + other.sigma**2)
return self
class GNSSStation(Location):
''' Representation of a GNSS station during a campaign measurement
For more information see
http://kb.unavco.org/kb/assets/660/UNAVCO_Campaign_GPS_GNSS_Handbook.pdf
'''
code = String.T(
help='Four letter station code',
optional=True)
style = StringChoice.T(
choices=['static', 'rapid_static', 'kinematic'],
default='static')
survey_start = DateTimestamp.T(
optional=True,
help='Survey start time')
survey_end = DateTimestamp.T(
optional=True,
help='Survey end time')
correlation_ne = Float.T(
default=0.,
help='North-East component correlation')
correlation_eu = Float.T(
default=0.,
help='East-Up component correlation')
correlation_nu = Float.T(
default=0.,
help='North-Up component correlation')
north = GNSSComponent.T(
optional=True)
east = GNSSComponent.T(
optional=True)
up = GNSSComponent.T(
optional=True)
def __eq__(self, other):
try:
return self.code == other.code
except AttributeError:
return False
def get_covariance_matrix(self):
components = self.components.values()
ncomponents = self.ncomponents
covar = num.zeros((ncomponents, ncomponents))
for ic1, comp1 in enumerate(components):
for ic2, comp2 in enumerate(components):
corr = self._get_comp_correlation(comp1, comp2)
covar[ic1, ic2] = corr * comp1.sigma * comp2.sigma
# This floating point operation is inaccurate:
# corr * comp1.sigma * comp2.sigma != corr * comp2.sigma * comp1.sigma
#
# Hence this identity
covar[num.tril_indices_from(covar, k=-1)] = \
covar[num.triu_indices_from(covar, k=1)]
return covar
def get_correlation_matrix(self):
components = self.components.values()
ncomponents = self.ncomponents
corr = num.zeros((ncomponents, ncomponents))
corr[num.diag_indices_from(corr)] = num.array(
[c.sigma for c in components])
for ic1, comp1 in enumerate(components):
for ic2, comp2 in enumerate(components):
if comp1 is comp2:
continue
corr[ic1, ic2] = self._get_comp_correlation(comp1, comp2)
# See comment at get_covariance_matrix
corr[num.tril_indices_from(corr, k=-1)] = \
corr[num.triu_indices_from(corr, k=1)]
return corr
def get_displacement_data(self):
return num.array([c.shift for c in self.components.values()])
def get_component_mask(self):
return num.array(
[False if self.__getattribute__(name) is None else True
for name in ('north', 'east', 'up')], dtype=num.bool)
@property
def components(self):
return OrderedDict(
[(name, self.__getattribute__(name))
for name in ('north', 'east', 'up')
if self.__getattribute__(name) is not None])
@property
def ncomponents(self):
return len(self.components)
def _get_comp_correlation(self, comp1, comp2):
if comp1 is comp2:
return 1.
s = self
correlation_map = {
(s.north, s.east): s.correlation_ne,
(s.east, s.up): s.correlation_eu,
(s.north, s.up): s.correlation_nu
}
return correlation_map.get(
(comp1, comp2),
correlation_map.get((comp2, comp1), False))
class GNSSCampaign(Object):
stations = List.T(
GNSSStation.T(),
help='List of GNSS campaign measurements')
name = String.T(
help='Campaign name',
default='Unnamed campaign')
survey_start = DateTimestamp.T(
optional=True)
survey_end = DateTimestamp.T(
optional=True)
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._cov_mat = None
self._cor_mat = None
def add_station(self, station):
self._cor_mat = None
self._cov_mat = None
return self.stations.append(station)
def remove_station(self, station_code):
try:
station = self.get_station(station_code)
self.stations.remove(station)
self._cor_mat = None
self._cov_mat = None
except ValueError:
logger.warn('Station {} does not exist in campaign, '
'do nothing.'.format(station_code))
def get_station(self, station_code):
for sta in self.stations:
if sta.code == station_code:
return sta
raise ValueError('Could not find station %s' % station_code)
def get_center_latlon(self):
return od.geographic_midpoint_locations(self.stations)
def get_radius(self):
coords = self.coordinates
return od.distance_accurate50m(
coords[:, 0].min(), coords[:, 1].min(),
coords[:, 0].max(), coords[:, 1].max()) / 2.
def get_covariance_matrix(self):
if self._cov_mat is None:
ncomponents = self.ncomponents
cov_arr = num.zeros((ncomponents, ncomponents))
idx = 0
for ista, sta in enumerate(self.stations):
ncomp = sta.ncomponents
cov_arr[idx:idx+ncomp, idx:idx+ncomp] = \
sta.get_covariance_matrix()
idx += ncomp
self._cov_mat = cov_arr
return self._cov_mat
def get_correlation_matrix(self):
if self._cor_mat is None:
ncomponents = self.ncomponents
cor_arr = num.zeros((ncomponents, ncomponents))
idx = 0
for ista, sta in enumerate(self.stations):
ncomp = sta.ncomponents
cor_arr[idx:idx+ncomp, idx:idx+ncomp] = \
sta.get_correlation_matrix()
idx += ncomp
self._cor_mat = cor_arr
return self._cor_mat
def get_component_mask(self):
return num.concatenate(
[s.get_component_mask() for s in self.stations])
def dump(self, *args, **kwargs):
self.regularize()
return Object.dump(self, *args, **kwargs)
@property
def coordinates(self):
return num.array([loc.effective_latlon for loc in self.stations])
@property
def nstations(self):
return len(self.stations)
@property
def ncomponents(self):
return sum([s.ncomponents for s in self.stations])
|