Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/gnss_campaign/target.py: 84%

1import logging 

2import numpy as num 

3from scipy import linalg as splinalg 

4 

5from pyrocko import gf 

6from pyrocko.guts import String, Dict, List, Int 

7from pyrocko.guts_array import Array 

8 

9from ..base import MisfitConfig, MisfitTarget, MisfitResult, TargetGroup 

10from grond.meta import has_get_plot_classes 

11 

12guts_prefix = 'grond' 

13logger = logging.getLogger('grond.target').getChild('gnss_campaign') 

14 

15 

16class GNSSCampaignMisfitResult(MisfitResult): 

17 """Carries the observations for a target and corresponding synthetics. """ 

18 statics_syn = Dict.T( 

19 String.T(), 

20 Array.T(dtype=num.float64, shape=(None,), serialize_as='list'), 

21 optional=True, 

22 help='Synthetic gnss surface displacements') 

23 statics_obs = Array.T( 

24 dtype=num.float64, shape=(None,), serialize_as='list', 

25 optional=True, 

26 help='Observed gnss surface displacements') 

27 

28 

29class GNSSCampaignMisfitConfig(MisfitConfig): 

30 pass 

31 

32 

33class GNSSCampaignTargetGroup(TargetGroup): 

34 """Handles static displacements from campaign GNSS observations, e.g GPS. 

35 

36 Station information, displacements and measurement errors are provided in 

37 a `yaml`-file (please find an example in the documentation). The 

38 measurement errors may consider correlations between components of a 

39 station, but correlations between single stations is not considered. 

40 """ 

41 gnss_campaigns = List.T( 

42 optional=True, 

43 help='List of individual campaign names' 

44 ' (`name` in `gnss.yaml` files).') 

45 misfit_config = GNSSCampaignMisfitConfig.T() 

46 

47 def get_targets(self, ds, event, default_path='none'): 

48 logger.debug('Selecting GNSS targets...') 

49 targets = [] 

50 

51 for camp in ds.get_gnss_campaigns(): 

52 if camp.name not in self.gnss_campaigns and\ 

53 '*all' not in self.gnss_campaigns: 

54 continue 

55 

56 if not isinstance(self.misfit_config, 

57 GNSSCampaignMisfitConfig): 

58 raise AttributeError('misfit_config must be of type' 

59 ' GNSSCampaignMisfitConfig.') 

60 

61 lats = num.array([s.lat for s in camp.stations]) 

62 lons = num.array([s.lon for s in camp.stations]) 

63 

64 north_shifts = num.array([s.north_shift for s in camp.stations]) 

65 east_shifts = num.array([s.east_shift for s in camp.stations]) 

66 

67 gnss_target = GNSSCampaignMisfitTarget( 

68 quantity='displacement', 

69 campaign_name=camp.name, 

70 lats=lats, 

71 lons=lons, 

72 east_shifts=east_shifts, 

73 north_shifts=north_shifts, 

74 ncomponents=camp.ncomponents, 

75 tsnapshot=None, 

76 interpolation=self.interpolation, 

77 store_id=self.store_id, 

78 normalisation_family=self.normalisation_family, 

79 path=self.path or default_path, 

80 misfit_config=self.misfit_config) 

81 

82 gnss_target.set_dataset(ds) 

83 targets.append(gnss_target) 

84 

85 return targets 

86 

87 

88@has_get_plot_classes 

89class GNSSCampaignMisfitTarget(gf.GNSSCampaignTarget, MisfitTarget): 

90 """Handles and carries out operations related to the objective functions. 

91 

92 The objective function is here the weighted misfit between observed 

93 and predicted surface displacements. 

94 """ 

95 campaign_name = String.T() 

96 ncomponents = Int.T(optional=True) 

97 misfit_config = GNSSCampaignMisfitConfig.T() 

98 

99 can_bootstrap_weights = True 

100 can_bootstrap_residuals = True 

101 

102 plot_misfits_cumulative = False 

103 

104 def __init__(self, **kwargs): 

105 gf.GNSSCampaignTarget.__init__(self, **kwargs) 

106 MisfitTarget.__init__(self, **kwargs) 

107 self._obs_data = None 

108 self._sigma = None 

109 self._weights = None 

110 self._correlated_weights = None 

111 self._station_component_mask = None 

112 

113 @property 

114 def id(self): 

115 return self.campaign_name 

116 

117 def string_id(self): 

118 return self.campaign_name 

119 

120 def misfits_string_ids(self): 

121 id_list = [] 

122 for station in self.campaign.stations: 

123 for name in ('north', 'east', 'up'): 

124 component = station.__getattribute__(name) 

125 if not component: 

126 continue 

127 id_list.append('%s.%s.%s' % 

128 (self.path, station.code, name[0].upper())) 

129 return id_list 

130 

131 @property 

132 def station_names(self): 

133 return ['%s' % (station.code) 

134 for station in self.campaign.stations] 

135 

136 @property 

137 def nmisfits(self): 

138 if self.ncomponents is None: 

139 try: 

140 self.campaign.ncomponents 

141 except ValueError: 

142 raise ValueError('Set the dataset!') 

143 return self.ncomponents 

144 

145 @property 

146 def nstations(self): 

147 return self.lats.size 

148 

149 def set_dataset(self, ds): 

150 MisfitTarget.set_dataset(self, ds) 

151 

152 def get_correlated_weights(self): 

153 if self._correlated_weights is None: 

154 self._correlated_weights = splinalg.sqrtm(self.weights) 

155 return self._correlated_weights 

156 

157 @property 

158 def campaign(self): 

159 return self._ds.get_gnss_campaign(self.campaign_name) 

160 

161 @property 

162 def obs_data(self): 

163 if self._obs_data is None: 

164 self._obs_data = num.concatenate( 

165 [s.get_displacement_data() for s in self.campaign.stations]) 

166 return self._obs_data 

167 

168 @property 

169 def obs_sigma(self): 

170 if self._sigma is None: 

171 self._sigma = num.array([ 

172 [s.north.sigma for s in self.campaign.stations], 

173 [s.east.sigma for s in self.campaign.stations], 

174 [s.up.sigma for s in self.campaign.stations]])\ 

175 .ravel(order='F') 

176 return self._sigma 

177 

178 @property 

179 def weights(self): 

180 """ 

181 Weights are the square-rooted, inverted data error variance-covariance. 

182 

183 The single component variances, and if provided the component 

184 covariances, are used to build a data variance matrix or 

185 variance-covariance matrix. 

186 

187 This matrix has the size for all possible NEU components, 

188 but throws zeros for not given components, also recorded in 

189 the _station_component_mask. 

190 """ 

191 if self._weights is None: 

192 covar = self.campaign.get_covariance_matrix() 

193 

194 if not num.any(covar.diagonal()): 

195 logger.warning('GNSS Stations have an empty covariance matrix.' 

196 ' Weights will be all equal.') 

197 num.fill_diagonal(covar, 1.) 

198 

199 self._weights = num.asmatrix(covar).I 

200 

201 return self._weights 

202 

203 @property 

204 def station_component_mask(self): 

205 if self._station_component_mask is None: 

206 self._station_component_mask = self.campaign.get_component_mask() 

207 return self._station_component_mask 

208 

209 @property 

210 def station_weights(self): 

211 weights = num.diag(self.weights) 

212 return num.mean([weights[0::3], weights[1::3], weights[2::3]], axis=0) 

213 

214 def component_weights(self): 

215 ws = num.sum(self.weights, axis=0) 

216 return ws 

217 

218 def post_process(self, engine, source, statics): 

219 """Applies the objective function. 

220 

221 As a result the weighted misfits are given and the observed and 

222 synthetic data. 

223 """ 

224 obs = self.obs_data 

225 

226 # All data is ordered in vectors as 

227 # S1_n, S1_e, S1_u, ..., Sn_n, Sn_e, Sn_u. Hence (.ravel(order='F')) 

228 syn = num.array([ 

229 statics['displacement.n'], 

230 statics['displacement.e'], 

231 -statics['displacement.d']])\ 

232 .ravel(order='F') 

233 

234 syn = syn[self.station_component_mask] 

235 res = obs - syn 

236 

237 misfit_value = res 

238 misfit_norm = obs 

239 

240 mf = num.vstack((misfit_value, misfit_norm)).T 

241 

242 result = GNSSCampaignMisfitResult( 

243 misfits=mf) 

244 if self._result_mode == 'full': 

245 result.statics_syn = statics 

246 result.statics_obs = obs 

247 

248 return result 

249 

250 def get_combined_weight(self): 

251 """A given manual weight in the configuration is applied.""" 

252 if self._combined_weight is None: 

253 self._combined_weight = num.full(self.nmisfits, self.manual_weight) 

254 

255 return self._combined_weight 

256 

257 def init_bootstrap_residuals(self, nbootstraps, rstate=None): 

258 logger.info('GNSS campaign %s, bootstrapping residuals' 

259 ' from measurement uncertainties ...' 

260 % self.campaign.name) 

261 if rstate is None: 

262 rstate = num.random.RandomState() 

263 

264 campaign = self.campaign 

265 bootstraps = num.empty((nbootstraps, self.ncomponents)) 

266 

267 sigmas = num.diag(num.sqrt(campaign.get_covariance_matrix())) 

268 if not num.all(sigmas): 

269 logger.warning('Bootstrapping GNSS stations is meaningless,' 

270 ' all station\'s sigma are 0.0!') 

271 

272 for ibs in range(nbootstraps): 

273 syn_noise = rstate.normal(scale=sigmas.ravel()) 

274 bootstraps[ibs, :] = syn_noise 

275 

276 self.set_bootstrap_residuals(bootstraps) 

277 

278 @classmethod 

279 def get_plot_classes(cls): 

280 from . import plot 

281 plots = super(GNSSCampaignMisfitTarget, cls).get_plot_classes() 

282 plots.extend(plot.get_plot_classes()) 

283 return plots