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

140 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-25 08:34 +0000

1import logging 

2import numpy as num 

3from scipy import linalg as splinalg 

4 

5from pyrocko import gf 

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

7from pyrocko.guts_array import Array 

8 

9import os 

10 

11from grond.meta import Parameter, has_get_plot_classes 

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

13 

14guts_prefix = 'grond' 

15logger = logging.getLogger('grond.targets.satellite.target') 

16 

17 

18class SatelliteMisfitConfig(MisfitConfig): 

19 """Carries the misfit configuration.""" 

20 optimise_orbital_ramp = Bool.T( 

21 default=True, 

22 help='Switch to account for a linear orbital ramp or not') 

23 ranges = Dict.T( 

24 String.T(), gf.Range.T(), 

25 default={'offset': '-0.5 .. 0.5', 

26 'ramp_east': '-1e-7 .. 1e-7', 

27 'ramp_north': '-1e-7 .. 1e-7' 

28 }, 

29 help='These parameters give bounds for an offset [m], a linear' 

30 ' gradient in east direction [m/m] and a linear gradient in north' 

31 ' direction [m/m]. Note, while the optimisation of these ramps' 

32 ' is individual for each target, the ranges set here are common' 

33 ' for all satellite targets.') 

34 

35 

36class SatelliteTargetGroup(TargetGroup): 

37 r""" 

38 Handles maps of static ground motion from satellite observations (InSAR) 

39 

40 The InSAR displacement maps post-processed by the `pyrocko` module `kite` 

41 are usually `Quadtree` downsampled (Jonsson, 2002). Each data point has a 

42 latitude, longitude, Line-of-sight displacement value [m] as well as an 

43 orientation and elevation angle, which define the Line-of-Sight. The data 

44 are associated with a weight matrix, which is the inverse of a full 

45 variance-covariance matrix (Sudhaus \& Jonsson, 2009). In principle, these 

46 data sets could stem from pixel offset maps. See also the documentation of 

47 the `kite` module. 

48 """ 

49 kite_scenes = List.T( 

50 optional=True, 

51 help='List of InSAR data files prepared \ 

52 by the ``pyrocko`` module ``kite``') 

53 misfit_config = SatelliteMisfitConfig.T( 

54 help='Settings for the objective function of these targets') 

55 

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

57 logger.debug('Selecting satellite targets...') 

58 targets = [] 

59 

60 for scene in ds.get_kite_scenes(): 

61 if scene.meta.scene_id not in self.kite_scenes and\ 

62 '*all' not in self.kite_scenes: 

63 continue 

64 

65 qt = scene.quadtree 

66 

67 lats = num.empty(qt.nleaves) 

68 lons = num.empty(qt.nleaves) 

69 lats.fill(qt.frame.llLat) 

70 lons.fill(qt.frame.llLon) 

71 

72 if qt.frame.isDegree(): 

73 logger.debug('Target "%s" is referenced in degree.' 

74 % scene.meta.scene_id) 

75 lons += qt.leaf_focal_points[:, 0] 

76 lats += qt.leaf_focal_points[:, 1] 

77 east_shifts = num.zeros_like(lats) 

78 north_shifts = num.zeros_like(lats) 

79 elif qt.frame.isMeter(): 

80 logger.debug('Target "%s" is referenced in meter.' 

81 % scene.meta.scene_id) 

82 east_shifts = qt.leaf_focal_points[:, 0] 

83 north_shifts = qt.leaf_focal_points[:, 1] 

84 else: 

85 assert False 

86 

87 sat_target = SatelliteMisfitTarget( 

88 quantity='displacement', 

89 scene_id=scene.meta.scene_id, 

90 lats=lats, 

91 lons=lons, 

92 east_shifts=east_shifts, 

93 north_shifts=north_shifts, 

94 theta=qt.leaf_thetas, 

95 phi=qt.leaf_phis, 

96 tsnapshot=None, 

97 interpolation=self.interpolation, 

98 store_id=self.store_id, 

99 normalisation_family=self.normalisation_family, 

100 path=self.path or default_path, 

101 misfit_config=self.misfit_config) 

102 

103 sat_target.set_dataset(ds) 

104 targets.append(sat_target) 

105 

106 return targets 

107 

108 

109class SatelliteMisfitResult(gf.Result, MisfitResult): 

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

111 statics_syn = Dict.T( 

112 String.T(), 

113 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'), 

114 optional=True, 

115 help='Predicted static displacements for a target (synthetics).') 

116 statics_obs = Dict.T( 

117 String.T(), 

118 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'), 

119 optional=True, 

120 help='Observed static displacement for a target.') 

121 

122 

123@has_get_plot_classes 

124class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget): 

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

126 

127 Standard operations are the calculation of the weighted misfit between 

128 observed and predicted (synthetic) data. If enabled in the misfit 

129 configuration, orbital ramps are optimized for. 

130 """ 

131 scene_id = String.T( 

132 help='UID string each Kite displacemente scene.' 

133 ' Corresponds to Kite scene_id.') 

134 misfit_config = SatelliteMisfitConfig.T( 

135 help='Configuration of the ``SatelliteTarget``') 

136 

137 can_bootstrap_residuals = True 

138 

139 available_parameters = [ 

140 Parameter('offset', 'm'), 

141 Parameter('ramp_north', 'm/m'), 

142 Parameter('ramp_east', 'm/m')] 

143 

144 def __init__(self, *args, **kwargs): 

145 gf.SatelliteTarget.__init__(self, *args, **kwargs) 

146 MisfitTarget.__init__(self, **kwargs) 

147 if not self.misfit_config.optimise_orbital_ramp: 

148 self.parameters = [] 

149 else: 

150 self.parameters = self.available_parameters 

151 

152 self.parameter_values = {} 

153 

154 self._noise_weight_matrix = None 

155 

156 @property 

157 def target_ranges(self): 

158 return self.misfit_config.ranges 

159 

160 def string_id(self): 

161 return '.'.join([self.path, self.scene_id]) 

162 

163 @property 

164 def id(self): 

165 return self.scene_id 

166 

167 def set_dataset(self, ds): 

168 MisfitTarget.set_dataset(self, ds) 

169 

170 @property 

171 def nmisfits(self): 

172 return self.lats.size 

173 

174 def get_correlated_weights(self, nthreads=0): 

175 ''' is for L2-norm weighting, the square-rooted, inverse covar ''' 

176 if self._noise_weight_matrix is None: 

177 logger.info( 

178 'Inverting scene covariance matrix (nthreads=%i)...' 

179 % nthreads) 

180 cov = self.scene.covariance 

181 cov.nthreads = nthreads 

182 

183 self._noise_weight_matrix = splinalg.sqrtm( 

184 num.linalg.inv(cov.covariance_matrix)) 

185 

186 logger.info('Inverting scene covariance matrix done.') 

187 

188 return self._noise_weight_matrix 

189 

190 @property 

191 def scene(self): 

192 return self._ds.get_kite_scene(self.scene_id) 

193 

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

195 """Applies the objective function. 

196 

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

198 synthetic data. For the satellite target the orbital ramp is 

199 calculated and applied here.""" 

200 scene = self.scene 

201 quadtree = scene.quadtree 

202 

203 obs = quadtree.leaf_medians 

204 

205 if self.misfit_config.optimise_orbital_ramp: 

206 stat_level = num.full_like(obs, self.parameter_values['offset']) 

207 

208 stat_level += (quadtree.leaf_center_distance[:, 0] 

209 * self.parameter_values['ramp_east']) 

210 stat_level += (quadtree.leaf_center_distance[:, 1] 

211 * self.parameter_values['ramp_north']) 

212 statics['displacement.los'] += stat_level 

213 

214 stat_syn = statics['displacement.los'] 

215 

216 res = obs - stat_syn 

217 

218 misfit_value = res 

219 misfit_norm = obs 

220 

221 mf = num.vstack([misfit_value, misfit_norm]).T 

222 result = SatelliteMisfitResult( 

223 misfits=mf) 

224 

225 if self._result_mode == 'full': 

226 result.statics_syn = statics 

227 result.statics_obs = {'displacement.los': obs} 

228 

229 return result 

230 

231 def get_combined_weight(self): 

232 if self._combined_weight is None: 

233 # invcov = self.scene.covariance.weight_matrix 

234 # self._combined_weight = invcov * self.manual_weight 

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

236 

237 return self._combined_weight 

238 

239 def prepare_modelling(self, engine, source, targets): 

240 return [self] 

241 

242 def finalize_modelling( 

243 self, engine, source, modelling_targets, modelling_results): 

244 return modelling_results[0] 

245 

246 def init_bootstrap_residuals(self, nbootstraps, rstate=None, nthreads=0): 

247 logger.info( 

248 'Scene "%s", initializing bootstrapping residuals from noise ' 

249 'pertubations...' % self.scene_id) 

250 

251 if rstate is None: 

252 rstate = num.random.RandomState() 

253 

254 scene = self.scene 

255 qt = scene.quadtree 

256 cov = scene.covariance 

257 bootstraps = num.zeros((nbootstraps, qt.nleaves)) 

258 

259 try: 

260 # TODO:mi Signal handler is not given back to the main task! 

261 # This is a python3.7 bug 

262 logger.warning('Using multi-threading for SatelliteTargets. ' 

263 'Python 3.7 is buggy and needs to be killed hard:' 

264 ' `killall grond`') 

265 from concurrent.futures import ThreadPoolExecutor 

266 nthreads = os.cpu_count() if not nthreads else nthreads 

267 

268 with ThreadPoolExecutor(max_workers=nthreads) as executor: 

269 res = executor.map( 

270 cov.getQuadtreeNoise, 

271 [rstate for _ in range(nbootstraps)]) 

272 

273 for ibs, bs in enumerate(res): 

274 bootstraps[ibs, :] = bs 

275 except ImportError: 

276 for ibs in range(nbootstraps): 

277 if not (ibs+1) % 5: 

278 logger.info('Calculating noise realisation %d/%d.' 

279 % (ibs+1, nbootstraps)) 

280 bootstraps[ibs, :] = cov.getQuadtreeNoise(rstate=rstate) 

281 

282 self.set_bootstrap_residuals(bootstraps) 

283 

284 @classmethod 

285 def get_plot_classes(cls): 

286 from . import plot 

287 plots = super(SatelliteMisfitTarget, cls).get_plot_classes() 

288 plots.extend(plot.get_plot_classes()) 

289 return plots 

290 

291 

292__all__ = ''' 

293 SatelliteTargetGroup 

294 SatelliteMisfitConfig 

295 SatelliteMisfitTarget 

296 SatelliteMisfitResult 

297'''.split()