Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/gnss_campaign/target.py: 84%
153 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:23 +0000
1import logging
2import numpy as num
3from scipy import linalg as splinalg
5from pyrocko import gf
6from pyrocko.guts import String, Dict, List, Int
8from ..base import MisfitConfig, MisfitTarget, MisfitResult, TargetGroup
9from grond.meta import has_get_plot_classes
11guts_prefix = 'grond'
12logger = logging.getLogger('grond.target').getChild('gnss_campaign')
15class GNSSCampaignMisfitResult(MisfitResult):
16 """Carries the observations for a target and corresponding synthetics. """
17 statics_syn = Dict.T(optional=True,
18 help='Synthetic gnss surface displacements')
19 statics_obs = Dict.T(optional=True,
20 help='Observed gnss surface displacements')
23class GNSSCampaignMisfitConfig(MisfitConfig):
24 pass
27class GNSSCampaignTargetGroup(TargetGroup):
28 """Handles static displacements from campaign GNSS observations, e.g GPS.
30 Station information, displacements and measurement errors are provided in
31 a `yaml`-file (please find an example in the documentation). The
32 measurement errors may consider correlations between components of a
33 station, but correlations between single stations is not considered.
34 """
35 gnss_campaigns = List.T(
36 optional=True,
37 help='List of individual campaign names'
38 ' (`name` in `gnss.yaml` files).')
39 misfit_config = GNSSCampaignMisfitConfig.T()
41 def get_targets(self, ds, event, default_path='none'):
42 logger.debug('Selecting GNSS targets...')
43 targets = []
45 for camp in ds.get_gnss_campaigns():
46 if camp.name not in self.gnss_campaigns and\
47 '*all' not in self.gnss_campaigns:
48 continue
50 if not isinstance(self.misfit_config,
51 GNSSCampaignMisfitConfig):
52 raise AttributeError('misfit_config must be of type'
53 ' GNSSCampaignMisfitConfig.')
55 lats = num.array([s.lat for s in camp.stations])
56 lons = num.array([s.lon for s in camp.stations])
58 north_shifts = num.array([s.north_shift for s in camp.stations])
59 east_shifts = num.array([s.east_shift for s in camp.stations])
61 gnss_target = GNSSCampaignMisfitTarget(
62 quantity='displacement',
63 campaign_name=camp.name,
64 lats=lats,
65 lons=lons,
66 east_shifts=east_shifts,
67 north_shifts=north_shifts,
68 ncomponents=camp.ncomponents,
69 tsnapshot=None,
70 interpolation=self.interpolation,
71 store_id=self.store_id,
72 normalisation_family=self.normalisation_family,
73 path=self.path or default_path,
74 misfit_config=self.misfit_config)
76 gnss_target.set_dataset(ds)
77 targets.append(gnss_target)
79 return targets
82@has_get_plot_classes
83class GNSSCampaignMisfitTarget(gf.GNSSCampaignTarget, MisfitTarget):
84 """Handles and carries out operations related to the objective functions.
86 The objective function is here the weighted misfit between observed
87 and predicted surface displacements.
88 """
89 campaign_name = String.T()
90 ncomponents = Int.T(optional=True)
91 misfit_config = GNSSCampaignMisfitConfig.T()
93 can_bootstrap_weights = True
94 can_bootstrap_residuals = True
96 plot_misfits_cumulative = False
98 def __init__(self, **kwargs):
99 gf.GNSSCampaignTarget.__init__(self, **kwargs)
100 MisfitTarget.__init__(self, **kwargs)
101 self._obs_data = None
102 self._sigma = None
103 self._weights = None
104 self._correlated_weights = None
105 self._station_component_mask = None
107 @property
108 def id(self):
109 return self.campaign_name
111 def string_id(self):
112 return self.campaign_name
114 def misfits_string_ids(self):
115 id_list = []
116 for station in self.campaign.stations:
117 for name in ('north', 'east', 'up'):
118 component = station.__getattribute__(name)
119 if not component:
120 continue
121 id_list.append('%s.%s.%s' %
122 (self.path, station.code, name[0].upper()))
123 return id_list
125 @property
126 def station_names(self):
127 return ['%s' % (station.code)
128 for station in self.campaign.stations]
130 @property
131 def nmisfits(self):
132 if self.ncomponents is None:
133 try:
134 self.campaign.ncomponents
135 except ValueError:
136 raise ValueError('Set the dataset!')
137 return self.ncomponents
139 @property
140 def nstations(self):
141 return self.lats.size
143 def set_dataset(self, ds):
144 MisfitTarget.set_dataset(self, ds)
146 def get_correlated_weights(self):
147 if self._correlated_weights is None:
148 self._correlated_weights = splinalg.sqrtm(self.weights)
149 return self._correlated_weights
151 @property
152 def campaign(self):
153 return self._ds.get_gnss_campaign(self.campaign_name)
155 @property
156 def obs_data(self):
157 if self._obs_data is None:
158 self._obs_data = num.concatenate(
159 [s.get_displacement_data() for s in self.campaign.stations])
160 return self._obs_data
162 @property
163 def obs_sigma(self):
164 if self._sigma is None:
165 self._sigma = num.array([
166 [s.north.sigma for s in self.campaign.stations],
167 [s.east.sigma for s in self.campaign.stations],
168 [s.up.sigma for s in self.campaign.stations]])\
169 .ravel(order='F')
170 return self._sigma
172 @property
173 def weights(self):
174 """
175 Weights are the square-rooted, inverted data error variance-covariance.
177 The single component variances, and if provided the component
178 covariances, are used to build a data variance matrix or
179 variance-covariance matrix.
181 This matrix has the size for all possible NEU components,
182 but throws zeros for not given components, also recorded in
183 the _station_component_mask.
184 """
185 if self._weights is None:
186 covar = self.campaign.get_covariance_matrix()
188 if not num.any(covar.diagonal()):
189 logger.warning('GNSS Stations have an empty covariance matrix.'
190 ' Weights will be all equal.')
191 num.fill_diagonal(covar, 1.)
193 self._weights = num.asmatrix(covar).I
195 return self._weights
197 @property
198 def station_component_mask(self):
199 if self._station_component_mask is None:
200 self._station_component_mask = self.campaign.get_component_mask()
201 return self._station_component_mask
203 @property
204 def station_weights(self):
205 weights = num.diag(self.weights)
206 return num.mean([weights[0::3], weights[1::3], weights[2::3]], axis=0)
208 def component_weights(self):
209 ws = num.sum(self.weights, axis=0)
210 return ws
212 def post_process(self, engine, source, statics):
213 """Applies the objective function.
215 As a result the weighted misfits are given and the observed and
216 synthetic data.
217 """
218 obs = self.obs_data
220 # All data is ordered in vectors as
221 # S1_n, S1_e, S1_u, ..., Sn_n, Sn_e, Sn_u. Hence (.ravel(order='F'))
222 syn = num.array([
223 statics['displacement.n'],
224 statics['displacement.e'],
225 -statics['displacement.d']])\
226 .ravel(order='F')
228 syn = syn[self.station_component_mask]
229 res = obs - syn
231 misfit_value = res
232 misfit_norm = obs
234 mf = num.vstack((misfit_value, misfit_norm)).T
236 result = GNSSCampaignMisfitResult(
237 misfits=mf)
238 if self._result_mode == 'full':
239 result.statics_syn = statics
240 result.statics_obs = obs
242 return result
244 def get_combined_weight(self):
245 """A given manual weight in the configuration is applied."""
246 if self._combined_weight is None:
247 self._combined_weight = num.full(self.nmisfits, self.manual_weight)
249 return self._combined_weight
251 def init_bootstrap_residuals(self, nbootstraps, rstate=None):
252 logger.info('GNSS campaign %s, bootstrapping residuals'
253 ' from measurement uncertainties ...'
254 % self.campaign.name)
255 if rstate is None:
256 rstate = num.random.RandomState()
258 campaign = self.campaign
259 bootstraps = num.empty((nbootstraps, self.ncomponents))
261 sigmas = num.diag(num.sqrt(campaign.get_covariance_matrix()))
262 if not num.all(sigmas):
263 logger.warning('Bootstrapping GNSS stations is meaningless,'
264 ' all station\'s sigma are 0.0!')
266 for ibs in range(nbootstraps):
267 syn_noise = rstate.normal(scale=sigmas.ravel())
268 bootstraps[ibs, :] = syn_noise
270 self.set_bootstrap_residuals(bootstraps)
272 @classmethod
273 def get_plot_classes(cls):
274 from . import plot
275 plots = super(GNSSCampaignMisfitTarget, cls).get_plot_classes()
276 plots.extend(plot.get_plot_classes())
277 return plots