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
6from pyrocko import gf
7from pyrocko.guts import String, Bool, Dict, List
8from pyrocko.guts_array import Array
11from grond.meta import Parameter, has_get_plot_classes
12from ..base import MisfitConfig, MisfitTarget, MisfitResult, TargetGroup
14guts_prefix = 'grond'
15logger = logging.getLogger('grond.targets.satellite.target')
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.')
36class SatelliteTargetGroup(TargetGroup):
37 r"""
38 Handles maps of static ground motion from satellite observations (InSAR)
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')
56 def get_targets(self, ds, event, default_path='none'):
57 logger.debug('Selecting satellite targets...')
58 targets = []
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
65 qt = scene.quadtree
67 lats = num.empty(qt.nleaves)
68 lons = num.empty(qt.nleaves)
69 lats.fill(qt.frame.llLat)
70 lons.fill(qt.frame.llLon)
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
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)
103 sat_target.set_dataset(ds)
104 targets.append(sat_target)
106 return targets
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.')
123@has_get_plot_classes
124class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
125 """Handles and carries out operations related to the objective functions.
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``')
137 can_bootstrap_residuals = True
139 available_parameters = [
140 Parameter('offset', 'm'),
141 Parameter('ramp_north', 'm/m'),
142 Parameter('ramp_east', 'm/m')]
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
152 self.parameter_values = {}
154 self._noise_weight_matrix = None
156 @property
157 def target_ranges(self):
158 return self.misfit_config.ranges
160 def string_id(self):
161 return '.'.join([self.path, self.scene_id])
163 @property
164 def id(self):
165 return self.scene_id
167 def set_dataset(self, ds):
168 MisfitTarget.set_dataset(self, ds)
170 @property
171 def nmisfits(self):
172 return self.lats.size
174 def get_correlated_weights(self):
175 ''' is for L2-norm weighting, the square-rooted, inverse covar '''
177 from grond.config import get_global_config
178 nthreads = get_global_config().nthreads
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
187 self._noise_weight_matrix = splinalg.sqrtm(
188 num.linalg.inv(cov.covariance_matrix))
190 logger.info('Inverting scene covariance matrix done.')
192 return self._noise_weight_matrix
194 @property
195 def scene(self):
196 return self._ds.get_kite_scene(self.scene_id)
198 def post_process(self, engine, source, statics):
199 """Applies the objective function.
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
207 obs = quadtree.leaf_medians
209 if self.misfit_config.optimise_orbital_ramp:
210 stat_level = num.full_like(obs, self.parameter_values['offset'])
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
218 stat_syn = statics['displacement.los']
220 res = obs - stat_syn
222 misfit_value = res
223 misfit_norm = obs
225 mf = num.vstack([misfit_value, misfit_norm]).T
226 result = SatelliteMisfitResult(
227 misfits=mf)
229 if self._result_mode == 'full':
230 result.statics_syn = statics
231 result.statics_obs = {'displacement.los': obs}
233 return result
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)
241 return self._combined_weight
243 def prepare_modelling(self, engine, source, targets):
244 return [self]
246 def finalize_modelling(
247 self, engine, source, modelling_targets, modelling_results):
248 return modelling_results[0]
250 def init_bootstrap_residuals(self, nbootstraps, rstate=None):
252 from grond.config import get_global_config
253 nthreads = get_global_config().nthreads
255 logger.info(
256 'Scene "%s", initializing bootstrapping residuals from noise '
257 'pertubations...' % self.scene_id)
259 if rstate is None:
260 rstate = num.random.RandomState()
262 scene = self.scene
263 qt = scene.quadtree
264 cov = scene.covariance
265 bootstraps = num.zeros((nbootstraps, qt.nleaves))
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`')
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)])
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)
291 self.set_bootstrap_residuals(bootstraps)
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
301__all__ = '''
302 SatelliteTargetGroup
303 SatelliteMisfitConfig
304 SatelliteMisfitTarget
305 SatelliteMisfitResult
306'''.split()