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

144 statements  

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

1import sys 

2import logging 

3import numpy as num 

4from scipy import linalg as splinalg 

5 

6from pyrocko import gf 

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

8from pyrocko.guts_array import Array 

9 

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

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

176 

177 from grond.config import get_global_config 

178 nthreads = get_global_config().nthreads 

179 

180 if self._noise_weight_matrix is None: 

181 logger.info( 

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

183 % nthreads) 

184 cov = self.scene.covariance 

185 cov.nthreads = nthreads 

186 

187 self._noise_weight_matrix = splinalg.sqrtm( 

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

189 

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

191 

192 return self._noise_weight_matrix 

193 

194 @property 

195 def scene(self): 

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

197 

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

199 """Applies the objective function. 

200 

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

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

203 calculated and applied here.""" 

204 scene = self.scene 

205 quadtree = scene.quadtree 

206 

207 obs = quadtree.leaf_medians 

208 

209 if self.misfit_config.optimise_orbital_ramp: 

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

211 

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

213 * self.parameter_values['ramp_east']) 

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

215 * self.parameter_values['ramp_north']) 

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

217 

218 stat_syn = statics['displacement.los'] 

219 

220 res = obs - stat_syn 

221 

222 misfit_value = res 

223 misfit_norm = obs 

224 

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

226 result = SatelliteMisfitResult( 

227 misfits=mf) 

228 

229 if self._result_mode == 'full': 

230 result.statics_syn = statics 

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

232 

233 return result 

234 

235 def get_combined_weight(self): 

236 if self._combined_weight is None: 

237 # invcov = self.scene.covariance.weight_matrix 

238 # self._combined_weight = invcov * self.manual_weight 

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

240 

241 return self._combined_weight 

242 

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

244 return [self] 

245 

246 def finalize_modelling( 

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

248 return modelling_results[0] 

249 

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

251 

252 from grond.config import get_global_config 

253 nthreads = get_global_config().nthreads 

254 

255 logger.info( 

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

257 'pertubations...' % self.scene_id) 

258 

259 if rstate is None: 

260 rstate = num.random.RandomState() 

261 

262 scene = self.scene 

263 qt = scene.quadtree 

264 cov = scene.covariance 

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

266 

267 try: 

268 if sys.version_info[:2] == (3, 7): 

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

270 # This is a python3.7 bug 

271 logger.warning( 

272 'Using multi-threading for SatelliteTargets. ' 

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

274 ' `killall grond`') 

275 

276 from concurrent.futures import ThreadPoolExecutor 

277 with ThreadPoolExecutor(max_workers=nthreads) as executor: 

278 res = executor.map( 

279 cov.getQuadtreeNoise, 

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

281 

282 for ibs, bs in enumerate(res): 

283 bootstraps[ibs, :] = bs 

284 except ImportError: 

285 for ibs in range(nbootstraps): 

286 if not (ibs+1) % 5: 

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

288 % (ibs+1, nbootstraps)) 

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

290 

291 self.set_bootstrap_residuals(bootstraps) 

292 

293 @classmethod 

294 def get_plot_classes(cls): 

295 from . import plot 

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

297 plots.extend(plot.get_plot_classes()) 

298 return plots 

299 

300 

301__all__ = ''' 

302 SatelliteTargetGroup 

303 SatelliteMisfitConfig 

304 SatelliteMisfitTarget 

305 SatelliteMisfitResult 

306'''.split()