Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/satellite/target.py: 94%
142 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 15:57 +0000
1import logging
2import numpy as num
3from scipy import linalg as splinalg
5from pyrocko import gf
6from pyrocko.guts import String, Bool, Dict, List
7from pyrocko.guts_array import Array
10from grond.meta import Parameter, has_get_plot_classes
11from ..base import MisfitConfig, MisfitTarget, MisfitResult, TargetGroup
13guts_prefix = 'grond'
14logger = logging.getLogger('grond.targets.satellite.target')
17class SatelliteMisfitConfig(MisfitConfig):
18 """Carries the misfit configuration."""
19 optimise_orbital_ramp = Bool.T(
20 default=True,
21 help='Switch to account for a linear orbital ramp or not')
22 ranges = Dict.T(
23 String.T(), gf.Range.T(),
24 default={'offset': '-0.5 .. 0.5',
25 'ramp_east': '-1e-7 .. 1e-7',
26 'ramp_north': '-1e-7 .. 1e-7'
27 },
28 help='These parameters give bounds for an offset [m], a linear'
29 ' gradient in east direction [m/m] and a linear gradient in north'
30 ' direction [m/m]. Note, while the optimisation of these ramps'
31 ' is individual for each target, the ranges set here are common'
32 ' for all satellite targets.')
35class SatelliteTargetGroup(TargetGroup):
36 r"""
37 Handles maps of static ground motion from satellite observations (InSAR)
39 The InSAR displacement maps post-processed by the `pyrocko` module `kite`
40 are usually `Quadtree` downsampled (Jonsson, 2002). Each data point has a
41 latitude, longitude, Line-of-sight displacement value [m] as well as an
42 orientation and elevation angle, which define the Line-of-Sight. The data
43 are associated with a weight matrix, which is the inverse of a full
44 variance-covariance matrix (Sudhaus \& Jonsson, 2009). In principle, these
45 data sets could stem from pixel offset maps. See also the documentation of
46 the `kite` module.
47 """
48 kite_scenes = List.T(
49 optional=True,
50 help='List of InSAR data files prepared \
51 by the ``pyrocko`` module ``kite``')
52 misfit_config = SatelliteMisfitConfig.T(
53 help='Settings for the objective function of these targets')
55 def get_targets(self, ds, event, default_path='none'):
56 logger.debug('Selecting satellite targets...')
57 targets = []
59 for scene in ds.get_kite_scenes():
60 if scene.meta.scene_id not in self.kite_scenes and\
61 '*all' not in self.kite_scenes:
62 continue
64 qt = scene.quadtree
66 lats = num.empty(qt.nleaves)
67 lons = num.empty(qt.nleaves)
68 lats.fill(qt.frame.llLat)
69 lons.fill(qt.frame.llLon)
71 if qt.frame.isDegree():
72 logger.debug('Target "%s" is referenced in degree.'
73 % scene.meta.scene_id)
74 lons += qt.leaf_focal_points[:, 0]
75 lats += qt.leaf_focal_points[:, 1]
76 east_shifts = num.zeros_like(lats)
77 north_shifts = num.zeros_like(lats)
78 elif qt.frame.isMeter():
79 logger.debug('Target "%s" is referenced in meter.'
80 % scene.meta.scene_id)
81 east_shifts = qt.leaf_focal_points[:, 0]
82 north_shifts = qt.leaf_focal_points[:, 1]
83 else:
84 assert False
86 sat_target = SatelliteMisfitTarget(
87 quantity='displacement',
88 scene_id=scene.meta.scene_id,
89 lats=lats,
90 lons=lons,
91 east_shifts=east_shifts,
92 north_shifts=north_shifts,
93 theta=qt.leaf_thetas,
94 phi=qt.leaf_phis,
95 tsnapshot=None,
96 interpolation=self.interpolation,
97 store_id=self.store_id,
98 normalisation_family=self.normalisation_family,
99 path=self.path or default_path,
100 misfit_config=self.misfit_config)
102 sat_target.set_dataset(ds)
103 targets.append(sat_target)
105 return targets
108class SatelliteMisfitResult(gf.Result, MisfitResult):
109 """Carries the observations for a target and corresponding synthetics."""
110 statics_syn = Dict.T(
111 String.T(),
112 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
113 optional=True,
114 help='Predicted static displacements for a target (synthetics).')
115 statics_obs = Dict.T(
116 String.T(),
117 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
118 optional=True,
119 help='Observed static displacement for a target.')
122@has_get_plot_classes
123class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
124 """Handles and carries out operations related to the objective functions.
126 Standard operations are the calculation of the weighted misfit between
127 observed and predicted (synthetic) data. If enabled in the misfit
128 configuration, orbital ramps are optimized for.
129 """
130 scene_id = String.T(
131 help='UID string each Kite displacemente scene.'
132 ' Corresponds to Kite scene_id.')
133 misfit_config = SatelliteMisfitConfig.T(
134 help='Configuration of the ``SatelliteTarget``')
136 can_bootstrap_residuals = True
138 available_parameters = [
139 Parameter('offset', 'm'),
140 Parameter('ramp_north', 'm/m'),
141 Parameter('ramp_east', 'm/m')]
143 def __init__(self, *args, **kwargs):
144 gf.SatelliteTarget.__init__(self, *args, **kwargs)
145 MisfitTarget.__init__(self, **kwargs)
146 if not self.misfit_config.optimise_orbital_ramp:
147 self.parameters = []
148 else:
149 self.parameters = self.available_parameters
151 self.parameter_values = {}
153 self._noise_weight_matrix = None
155 @property
156 def target_ranges(self):
157 return self.misfit_config.ranges
159 def string_id(self):
160 return '.'.join([self.path, self.scene_id])
162 @property
163 def id(self):
164 return self.scene_id
166 def set_dataset(self, ds):
167 MisfitTarget.set_dataset(self, ds)
169 @property
170 def nmisfits(self):
171 return self.lats.size
173 def get_correlated_weights(self):
174 ''' is for L2-norm weighting, the square-rooted, inverse covar '''
176 from grond.config import get_global_config
177 nthreads = get_global_config().nthreads
179 if self._noise_weight_matrix is None:
180 logger.info(
181 'Inverting scene covariance matrix (nthreads=%i)...'
182 % nthreads)
183 cov = self.scene.covariance
184 cov.nthreads = nthreads
186 self._noise_weight_matrix = splinalg.sqrtm(
187 num.linalg.inv(cov.covariance_matrix))
189 logger.info('Inverting scene covariance matrix done.')
191 return self._noise_weight_matrix
193 @property
194 def scene(self):
195 return self._ds.get_kite_scene(self.scene_id)
197 def post_process(self, engine, source, statics):
198 """Applies the objective function.
200 As a result the weighted misfits are given and the observed and
201 synthetic data. For the satellite target the orbital ramp is
202 calculated and applied here."""
203 scene = self.scene
204 quadtree = scene.quadtree
206 obs = quadtree.leaf_medians
208 if self.misfit_config.optimise_orbital_ramp:
209 stat_level = num.full_like(obs, self.parameter_values['offset'])
211 stat_level += (quadtree.leaf_center_distance[:, 0]
212 * self.parameter_values['ramp_east'])
213 stat_level += (quadtree.leaf_center_distance[:, 1]
214 * self.parameter_values['ramp_north'])
215 statics['displacement.los'] += stat_level
217 stat_syn = statics['displacement.los']
219 res = obs - stat_syn
221 misfit_value = res
222 misfit_norm = obs
224 mf = num.vstack([misfit_value, misfit_norm]).T
225 result = SatelliteMisfitResult(
226 misfits=mf)
228 if self._result_mode == 'full':
229 result.statics_syn = statics
230 result.statics_obs = {'displacement.los': obs}
232 return result
234 def get_combined_weight(self):
235 if self._combined_weight is None:
236 # invcov = self.scene.covariance.weight_matrix
237 # self._combined_weight = invcov * self.manual_weight
238 self._combined_weight = num.full(self.nmisfits, self.manual_weight)
240 return self._combined_weight
242 def prepare_modelling(self, engine, source, targets):
243 return [self]
245 def finalize_modelling(
246 self, engine, source, modelling_targets, modelling_results):
247 return modelling_results[0]
249 def init_bootstrap_residuals(self, nbootstraps, rstate=None):
251 from grond.config import get_global_config
252 nthreads = get_global_config().nthreads
254 logger.info(
255 'Scene "%s", initializing bootstrapping residuals from noise '
256 'pertubations...' % self.scene_id)
258 if rstate is None:
259 rstate = num.random.RandomState()
261 scene = self.scene
262 qt = scene.quadtree
263 cov = scene.covariance
264 bootstraps = num.zeros((nbootstraps, qt.nleaves))
266 try:
267 # TODO:mi Signal handler is not given back to the main task!
268 # This is a python3.7 bug
269 logger.warning('Using multi-threading for SatelliteTargets. '
270 'Python 3.7 is buggy and needs to be killed hard:'
271 ' `killall grond`')
272 from concurrent.futures import ThreadPoolExecutor
273 with ThreadPoolExecutor(max_workers=nthreads) as executor:
274 res = executor.map(
275 cov.getQuadtreeNoise,
276 [rstate for _ in range(nbootstraps)])
278 for ibs, bs in enumerate(res):
279 bootstraps[ibs, :] = bs
280 except ImportError:
281 for ibs in range(nbootstraps):
282 if not (ibs+1) % 5:
283 logger.info('Calculating noise realisation %d/%d.'
284 % (ibs+1, nbootstraps))
285 bootstraps[ibs, :] = cov.getQuadtreeNoise(rstate=rstate)
287 self.set_bootstrap_residuals(bootstraps)
289 @classmethod
290 def get_plot_classes(cls):
291 from . import plot
292 plots = super(SatelliteMisfitTarget, cls).get_plot_classes()
293 plots.extend(plot.get_plot_classes())
294 return plots
297__all__ = '''
298 SatelliteTargetGroup
299 SatelliteMisfitConfig
300 SatelliteMisfitTarget
301 SatelliteMisfitResult
302'''.split()