Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/satellite/target.py: 92%
145 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +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.')
33 use_cross_correlation = Bool.T(
34 default=False,
35 help='If ``True``, only the CC coefficient between observed and '
36 'modelled displacement is taken into account.')
39class SatelliteTargetGroup(TargetGroup):
40 r"""
41 Handles maps of static ground motion from satellite observations (InSAR)
43 The InSAR displacement maps post-processed by the `pyrocko` module `kite`
44 are usually `Quadtree` downsampled (Jonsson, 2002). Each data point has a
45 latitude, longitude, Line-of-sight displacement value [m] as well as an
46 orientation and elevation angle, which define the Line-of-Sight. The data
47 are associated with a weight matrix, which is the inverse of a full
48 variance-covariance matrix (Sudhaus \& Jonsson, 2009). In principle, these
49 data sets could stem from pixel offset maps. See also the documentation of
50 the `kite` module.
51 """
52 kite_scenes = List.T(
53 optional=True,
54 help='List of InSAR data files prepared \
55 by the ``pyrocko`` module ``kite``')
56 misfit_config = SatelliteMisfitConfig.T(
57 help='Settings for the objective function of these targets')
59 def get_targets(self, ds, event, default_path='none'):
60 logger.debug('Selecting satellite targets...')
61 targets = []
63 for scene in ds.get_kite_scenes():
64 if scene.meta.scene_id not in self.kite_scenes and\
65 '*all' not in self.kite_scenes:
66 continue
68 qt = scene.quadtree
70 lats = num.empty(qt.nleaves)
71 lons = num.empty(qt.nleaves)
72 lats.fill(qt.frame.llLat)
73 lons.fill(qt.frame.llLon)
75 if qt.frame.isDegree():
76 logger.debug('Target "%s" is referenced in degree.'
77 % scene.meta.scene_id)
78 lons += qt.leaf_focal_points[:, 0]
79 lats += qt.leaf_focal_points[:, 1]
80 east_shifts = num.zeros_like(lats)
81 north_shifts = num.zeros_like(lats)
82 elif qt.frame.isMeter():
83 logger.debug('Target "%s" is referenced in meter.'
84 % scene.meta.scene_id)
85 east_shifts = qt.leaf_focal_points[:, 0]
86 north_shifts = qt.leaf_focal_points[:, 1]
87 else:
88 assert False
90 sat_target = SatelliteMisfitTarget(
91 quantity='displacement',
92 scene_id=scene.meta.scene_id,
93 lats=lats,
94 lons=lons,
95 east_shifts=east_shifts,
96 north_shifts=north_shifts,
97 theta=qt.leaf_thetas,
98 phi=qt.leaf_phis,
99 tsnapshot=None,
100 interpolation=self.interpolation,
101 store_id=self.store_id,
102 normalisation_family=self.normalisation_family,
103 path=self.path or default_path,
104 misfit_config=self.misfit_config)
106 sat_target.set_dataset(ds)
107 targets.append(sat_target)
109 return targets
112class SatelliteMisfitResult(gf.Result, MisfitResult):
113 """Carries the observations for a target and corresponding synthetics."""
114 statics_syn = Dict.T(
115 String.T(),
116 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
117 optional=True,
118 help='Predicted static displacements for a target (synthetics).')
119 statics_obs = Dict.T(
120 String.T(),
121 Array.T(dtype=num.float64, shape=(None,), serialize_as='base64'),
122 optional=True,
123 help='Observed static displacement for a target.')
126@has_get_plot_classes
127class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
128 """Handles and carries out operations related to the objective functions.
130 Standard operations are the calculation of the weighted misfit between
131 observed and predicted (synthetic) data. If enabled in the misfit
132 configuration, orbital ramps are optimized for.
133 """
134 scene_id = String.T(
135 help='UID string each Kite displacemente scene.'
136 ' Corresponds to Kite scene_id.')
137 misfit_config = SatelliteMisfitConfig.T(
138 help='Configuration of the ``SatelliteTarget``')
140 can_bootstrap_residuals = True
142 available_parameters = [
143 Parameter('offset', 'm'),
144 Parameter('ramp_north', 'm/m'),
145 Parameter('ramp_east', 'm/m')]
147 def __init__(self, *args, **kwargs):
148 gf.SatelliteTarget.__init__(self, *args, **kwargs)
149 MisfitTarget.__init__(self, **kwargs)
150 if not self.misfit_config.optimise_orbital_ramp:
151 self.parameters = []
152 else:
153 self.parameters = self.available_parameters
155 self.parameter_values = {}
157 self._noise_weight_matrix = None
159 @property
160 def target_ranges(self):
161 return self.misfit_config.ranges
163 def string_id(self):
164 return '.'.join([self.path, self.scene_id])
166 @property
167 def id(self):
168 return self.scene_id
170 def set_dataset(self, ds):
171 MisfitTarget.set_dataset(self, ds)
173 @property
174 def nmisfits(self):
175 return self.lats.size
177 def get_correlated_weights(self):
178 ''' is for L2-norm weighting, the square-rooted, inverse covar '''
180 from grond.config import get_global_config
181 nthreads = get_global_config().nthreads
183 if self._noise_weight_matrix is None:
184 logger.info(
185 'Inverting scene covariance matrix (nthreads=%i)...'
186 % nthreads)
187 cov = self.scene.covariance
188 cov.nthreads = nthreads
190 self._noise_weight_matrix = splinalg.sqrtm(
191 num.linalg.inv(cov.covariance_matrix))
193 logger.info('Inverting scene covariance matrix done.')
195 return self._noise_weight_matrix
197 @property
198 def scene(self):
199 return self._ds.get_kite_scene(self.scene_id)
201 def post_process(self, engine, source, statics):
202 """Applies the objective function.
204 As a result the weighted misfits are given and the observed and
205 synthetic data. For the satellite target the orbital ramp is
206 calculated and applied here."""
207 scene = self.scene
208 quadtree = scene.quadtree
210 obs = quadtree.leaf_medians
212 if self.misfit_config.optimise_orbital_ramp:
213 stat_level = num.full_like(obs, self.parameter_values['offset'])
215 stat_level += (quadtree.leaf_center_distance[:, 0]
216 * self.parameter_values['ramp_east'])
217 stat_level += (quadtree.leaf_center_distance[:, 1]
218 * self.parameter_values['ramp_north'])
219 statics['displacement.los'] += stat_level
221 stat_syn = statics['displacement.los']
223 if self.misfit_config.use_cross_correlation:
224 obs /= num.abs(obs).max()
225 stat_syn /= num.abs(stat_syn).max()
227 res = obs - stat_syn
229 misfit_value = res
230 misfit_norm = obs
232 mf = num.vstack([misfit_value, misfit_norm]).T
233 result = SatelliteMisfitResult(
234 misfits=mf)
236 if self._result_mode == 'full':
237 result.statics_syn = statics
238 result.statics_obs = {'displacement.los': obs}
240 return result
242 def get_combined_weight(self):
243 if self._combined_weight is None:
244 # invcov = self.scene.covariance.weight_matrix
245 # self._combined_weight = invcov * self.manual_weight
246 self._combined_weight = num.full(self.nmisfits, self.manual_weight)
248 return self._combined_weight
250 def prepare_modelling(self, engine, source, targets):
251 return [self]
253 def finalize_modelling(
254 self, engine, source, modelling_targets, modelling_results):
255 return modelling_results[0]
257 def init_bootstrap_residuals(self, nbootstraps, rstate=None):
259 from grond.config import get_global_config
260 nthreads = get_global_config().nthreads
262 logger.info(
263 'Scene "%s", initializing bootstrapping residuals from noise '
264 'pertubations...' % self.scene_id)
266 if rstate is None:
267 rstate = num.random.RandomState()
269 scene = self.scene
270 qt = scene.quadtree
271 cov = scene.covariance
272 bootstraps = num.zeros((nbootstraps, qt.nleaves))
274 try:
275 # TODO:mi Signal handler is not given back to the main task!
276 # This is a python3.7 bug
277 from concurrent.futures import ThreadPoolExecutor
278 with ThreadPoolExecutor(max_workers=nthreads) as executor:
279 res = executor.map(
280 cov.getQuadtreeNoise,
281 [rstate for _ in range(nbootstraps)])
283 for ibs, bs in enumerate(res):
284 bootstraps[ibs, :] = bs
285 except ImportError:
286 for ibs in range(nbootstraps):
287 if not (ibs+1) % 5:
288 logger.info('Calculating noise realisation %d/%d.'
289 % (ibs+1, nbootstraps))
290 bootstraps[ibs, :] = cov.getQuadtreeNoise(rstate=rstate)
292 self.set_bootstrap_residuals(bootstraps)
294 @classmethod
295 def get_plot_classes(cls):
296 from . import plot
297 plots = super(SatelliteMisfitTarget, cls).get_plot_classes()
298 plots.extend(plot.get_plot_classes())
299 return plots
302__all__ = '''
303 SatelliteTargetGroup
304 SatelliteMisfitConfig
305 SatelliteMisfitTarget
306 SatelliteMisfitResult
307'''.split()