Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/satellite/target.py: 92%
144 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-25 10:12 +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
9import os
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.')
34 use_cross_correlation = Bool.T(
35 default=False,
36 help='If ``True``, only the CC coefficient between observed and '
37 'modelled displacement is taken into account.')
40class SatelliteTargetGroup(TargetGroup):
41 r"""
42 Handles maps of static ground motion from satellite observations (InSAR)
44 The InSAR displacement maps post-processed by the `pyrocko` module `kite`
45 are usually `Quadtree` downsampled (Jonsson, 2002). Each data point has a
46 latitude, longitude, Line-of-sight displacement value [m] as well as an
47 orientation and elevation angle, which define the Line-of-Sight. The data
48 are associated with a weight matrix, which is the inverse of a full
49 variance-covariance matrix (Sudhaus \& Jonsson, 2009). In principle, these
50 data sets could stem from pixel offset maps. See also the documentation of
51 the `kite` module.
52 """
53 kite_scenes = List.T(
54 optional=True,
55 help='List of InSAR data files prepared \
56 by the ``pyrocko`` module ``kite``')
57 misfit_config = SatelliteMisfitConfig.T(
58 help='Settings for the objective function of these targets')
60 def get_targets(self, ds, event, default_path='none'):
61 logger.debug('Selecting satellite targets...')
62 targets = []
64 for scene in ds.get_kite_scenes():
65 if scene.meta.scene_id not in self.kite_scenes and\
66 '*all' not in self.kite_scenes:
67 continue
69 qt = scene.quadtree
71 lats = num.empty(qt.nleaves)
72 lons = num.empty(qt.nleaves)
73 lats.fill(qt.frame.llLat)
74 lons.fill(qt.frame.llLon)
76 if qt.frame.isDegree():
77 logger.debug('Target "%s" is referenced in degree.'
78 % scene.meta.scene_id)
79 lons += qt.leaf_focal_points[:, 0]
80 lats += qt.leaf_focal_points[:, 1]
81 east_shifts = num.zeros_like(lats)
82 north_shifts = num.zeros_like(lats)
83 elif qt.frame.isMeter():
84 logger.debug('Target "%s" is referenced in meter.'
85 % scene.meta.scene_id)
86 east_shifts = qt.leaf_focal_points[:, 0]
87 north_shifts = qt.leaf_focal_points[:, 1]
88 else:
89 assert False
91 sat_target = SatelliteMisfitTarget(
92 quantity='displacement',
93 scene_id=scene.meta.scene_id,
94 lats=lats,
95 lons=lons,
96 east_shifts=east_shifts,
97 north_shifts=north_shifts,
98 theta=qt.leaf_thetas,
99 phi=qt.leaf_phis,
100 tsnapshot=None,
101 interpolation=self.interpolation,
102 store_id=self.store_id,
103 normalisation_family=self.normalisation_family,
104 path=self.path or default_path,
105 misfit_config=self.misfit_config)
107 sat_target.set_dataset(ds)
108 targets.append(sat_target)
110 return targets
113class SatelliteMisfitResult(gf.Result, MisfitResult):
114 """Carries the observations for a target and corresponding synthetics."""
115 statics_syn = Dict.T(
116 String.T(),
117 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
118 optional=True,
119 help='Predicted static displacements for a target (synthetics).')
120 statics_obs = Dict.T(
121 String.T(),
122 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
123 optional=True,
124 help='Observed static displacement for a target.')
127@has_get_plot_classes
128class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
129 """Handles and carries out operations related to the objective functions.
131 Standard operations are the calculation of the weighted misfit between
132 observed and predicted (synthetic) data. If enabled in the misfit
133 configuration, orbital ramps are optimized for.
134 """
135 scene_id = String.T(
136 help='UID string each Kite displacemente scene.'
137 ' Corresponds to Kite scene_id.')
138 misfit_config = SatelliteMisfitConfig.T(
139 help='Configuration of the ``SatelliteTarget``')
141 can_bootstrap_residuals = True
143 available_parameters = [
144 Parameter('offset', 'm'),
145 Parameter('ramp_north', 'm/m'),
146 Parameter('ramp_east', 'm/m')]
148 def __init__(self, *args, **kwargs):
149 gf.SatelliteTarget.__init__(self, *args, **kwargs)
150 MisfitTarget.__init__(self, **kwargs)
151 if not self.misfit_config.optimise_orbital_ramp:
152 self.parameters = []
153 else:
154 self.parameters = self.available_parameters
156 self.parameter_values = {}
158 self._noise_weight_matrix = None
160 @property
161 def target_ranges(self):
162 return self.misfit_config.ranges
164 def string_id(self):
165 return '.'.join([self.path, self.scene_id])
167 @property
168 def id(self):
169 return self.scene_id
171 def set_dataset(self, ds):
172 MisfitTarget.set_dataset(self, ds)
174 @property
175 def nmisfits(self):
176 return self.lats.size
178 def get_correlated_weights(self, nthreads=0):
179 ''' is for L2-norm weighting, the square-rooted, inverse covar '''
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 if self.misfit_config.use_cross_correlation:
221 obs /= num.abs(obs).max()
222 stat_syn /= num.abs(stat_syn).max()
224 res = obs - stat_syn
226 misfit_value = res
227 misfit_norm = obs
229 mf = num.vstack([misfit_value, misfit_norm]).T
230 result = SatelliteMisfitResult(
231 misfits=mf)
233 if self._result_mode == 'full':
234 result.statics_syn = statics
235 result.statics_obs = {'displacement.los': obs}
237 return result
239 def get_combined_weight(self):
240 if self._combined_weight is None:
241 # invcov = self.scene.covariance.weight_matrix
242 # self._combined_weight = invcov * self.manual_weight
243 self._combined_weight = num.full(self.nmisfits, self.manual_weight)
245 return self._combined_weight
247 def prepare_modelling(self, engine, source, targets):
248 return [self]
250 def finalize_modelling(
251 self, engine, source, modelling_targets, modelling_results):
252 return modelling_results[0]
254 def init_bootstrap_residuals(self, nbootstraps, rstate=None, nthreads=0):
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 # TODO:mi Signal handler is not given back to the main task!
269 # This is a python3.7 bug
270 from concurrent.futures import ThreadPoolExecutor
271 nthreads = os.cpu_count() if not nthreads else nthreads
272 logger.debug('Using %d threads'
273 ' for bootstrapping SatelliteTargets.', nthreads)
275 with ThreadPoolExecutor(max_workers=nthreads) as executor:
276 res = executor.map(
277 cov.getQuadtreeNoise,
278 [rstate for _ in range(nbootstraps)])
280 for ibs, bs in enumerate(res):
281 bootstraps[ibs, :] = bs
282 except ImportError:
283 for ibs in range(nbootstraps):
284 if not (ibs+1) % 5:
285 logger.info('Calculating noise realisation %d/%d.'
286 % (ibs+1, nbootstraps))
287 bootstraps[ibs, :] = cov.getQuadtreeNoise(rstate=rstate)
289 self.set_bootstrap_residuals(bootstraps)
291 @classmethod
292 def get_plot_classes(cls):
293 from . import plot
294 plots = super(SatelliteMisfitTarget, cls).get_plot_classes()
295 plots.extend(plot.get_plot_classes())
296 return plots
299__all__ = '''
300 SatelliteTargetGroup
301 SatelliteMisfitConfig
302 SatelliteMisfitTarget
303 SatelliteMisfitResult
304'''.split()