Coverage for /usr/local/lib/python3.7/dist-packages/pyrocko/model/gnss.py : 38%

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])