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

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import logging 

5import numpy as num 

6from scipy import linalg as splinalg 

7 

8from pyrocko import gf 

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

10 

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

12from grond.meta import has_get_plot_classes 

13 

14guts_prefix = 'grond' 

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

16 

17 

18class GNSSCampaignMisfitResult(MisfitResult): 

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

20 statics_syn = Dict.T(optional=True, 

21 help='Synthetic gnss surface displacements') 

22 statics_obs = Dict.T(optional=True, 

23 help='Observed gnss surface displacements') 

24 

25 

26class GNSSCampaignMisfitConfig(MisfitConfig): 

27 pass 

28 

29 

30class GNSSCampaignTargetGroup(TargetGroup): 

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

32 

33 Station information, displacements and measurement errors are provided in 

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

35 measurement errors may consider correlations between components of a 

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

37 """ 

38 gnss_campaigns = List.T( 

39 optional=True, 

40 help='List of individual campaign names' 

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

42 misfit_config = GNSSCampaignMisfitConfig.T() 

43 

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

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

46 targets = [] 

47 

48 for camp in ds.get_gnss_campaigns(): 

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

50 '*all' not in self.gnss_campaigns: 

51 continue 

52 

53 if not isinstance(self.misfit_config, 

54 GNSSCampaignMisfitConfig): 

55 raise AttributeError('misfit_config must be of type' 

56 ' GNSSCampaignMisfitConfig.') 

57 

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

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

60 

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

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

63 

64 gnss_target = GNSSCampaignMisfitTarget( 

65 quantity='displacement', 

66 campaign_name=camp.name, 

67 lats=lats, 

68 lons=lons, 

69 east_shifts=east_shifts, 

70 north_shifts=north_shifts, 

71 ncomponents=camp.ncomponents, 

72 tsnapshot=None, 

73 interpolation=self.interpolation, 

74 store_id=self.store_id, 

75 normalisation_family=self.normalisation_family, 

76 path=self.path or default_path, 

77 misfit_config=self.misfit_config) 

78 

79 gnss_target.set_dataset(ds) 

80 targets.append(gnss_target) 

81 

82 return targets 

83 

84 

85@has_get_plot_classes 

86class GNSSCampaignMisfitTarget(gf.GNSSCampaignTarget, MisfitTarget): 

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

88 

89 The objective function is here the weighted misfit between observed 

90 and predicted surface displacements. 

91 """ 

92 campaign_name = String.T() 

93 ncomponents = Int.T(optional=True) 

94 misfit_config = GNSSCampaignMisfitConfig.T() 

95 

96 can_bootstrap_weights = True 

97 can_bootstrap_residuals = True 

98 

99 plot_misfits_cumulative = False 

100 

101 def __init__(self, **kwargs): 

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

103 MisfitTarget.__init__(self, **kwargs) 

104 self._obs_data = None 

105 self._sigma = None 

106 self._weights = None 

107 self._correlated_weights = None 

108 self._station_component_mask = None 

109 

110 @property 

111 def id(self): 

112 return self.campaign_name 

113 

114 def string_id(self): 

115 return self.campaign_name 

116 

117 def misfits_string_ids(self): 

118 id_list = [] 

119 for station in self.campaign.stations: 

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

121 component = station.__getattribute__(name) 

122 if not component: 

123 continue 

124 id_list.append('%s.%s.%s' % 

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

126 return id_list 

127 

128 @property 

129 def station_names(self): 

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

131 for station in self.campaign.stations] 

132 

133 @property 

134 def nmisfits(self): 

135 if self.ncomponents is None: 

136 try: 

137 self.campaign.ncomponents 

138 except ValueError: 

139 raise ValueError('Set the dataset!') 

140 return self.ncomponents 

141 

142 @property 

143 def nstations(self): 

144 return self.lats.size 

145 

146 def set_dataset(self, ds): 

147 MisfitTarget.set_dataset(self, ds) 

148 

149 def get_correlated_weights(self): 

150 if self._correlated_weights is None: 

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

152 return self._correlated_weights 

153 

154 @property 

155 def campaign(self): 

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

157 

158 @property 

159 def obs_data(self): 

160 if self._obs_data is None: 

161 self._obs_data = num.concatenate( 

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

163 return self._obs_data 

164 

165 @property 

166 def obs_sigma(self): 

167 if self._sigma is None: 

168 self._sigma = num.array([ 

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

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

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

172 .ravel(order='F') 

173 return self._sigma 

174 

175 @property 

176 def weights(self): 

177 """ 

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

179 

180 The single component variances, and if provided the component 

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

182 variance-covariance matrix. 

183 

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

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

186 the _station_component_mask. 

187 """ 

188 if self._weights is None: 

189 covar = self.campaign.get_covariance_matrix() 

190 

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

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

193 ' Weights will be all equal.') 

194 num.fill_diagonal(covar, 1.) 

195 

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

197 

198 return self._weights 

199 

200 @property 

201 def station_component_mask(self): 

202 if self._station_component_mask is None: 

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

204 return self._station_component_mask 

205 

206 @property 

207 def station_weights(self): 

208 weights = num.diag(self.weights) 

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

210 

211 def component_weights(self): 

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

213 return ws 

214 

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

216 """Applies the objective function. 

217 

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

219 synthetic data. 

220 """ 

221 obs = self.obs_data 

222 

223 # All data is ordered in vectors as 

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

225 syn = num.array([ 

226 statics['displacement.n'], 

227 statics['displacement.e'], 

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

229 .ravel(order='F') 

230 

231 syn = syn[self.station_component_mask] 

232 res = obs - syn 

233 

234 misfit_value = res 

235 misfit_norm = obs 

236 

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

238 

239 result = GNSSCampaignMisfitResult( 

240 misfits=mf) 

241 if self._result_mode == 'full': 

242 result.statics_syn = statics 

243 result.statics_obs = obs 

244 

245 return result 

246 

247 def get_combined_weight(self): 

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

249 if self._combined_weight is None: 

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

251 

252 return self._combined_weight 

253 

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

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

256 ' from measurement uncertainties ...' 

257 % self.campaign.name) 

258 if rstate is None: 

259 rstate = num.random.RandomState() 

260 

261 campaign = self.campaign 

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

263 

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

265 if not num.all(sigmas): 

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

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

268 

269 for ibs in range(nbootstraps): 

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

271 bootstraps[ibs, :] = syn_noise 

272 

273 self.set_bootstrap_residuals(bootstraps) 

274 

275 @classmethod 

276 def get_plot_classes(cls): 

277 from . import plot 

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

279 plots.extend(plot.get_plot_classes()) 

280 return plots