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.')
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, nthreads=0):
175 ''' is for L2-norm weighting, the square-rooted, inverse covar '''
176 if self._noise_weight_matrix is None:
177 logger.info(
178 'Inverting scene covariance matrix (nthreads=%i)...'
179 % nthreads)
180 cov = self.scene.covariance
181 cov.nthreads = nthreads
183 self._noise_weight_matrix = splinalg.sqrtm(
184 num.linalg.inv(cov.covariance_matrix))
186 logger.info('Inverting scene covariance matrix done.')
188 return self._noise_weight_matrix
190 @property
191 def scene(self):
192 return self._ds.get_kite_scene(self.scene_id)
194 def post_process(self, engine, source, statics):
195 """Applies the objective function.
197 As a result the weighted misfits are given and the observed and
198 synthetic data. For the satellite target the orbital ramp is
199 calculated and applied here."""
200 scene = self.scene
201 quadtree = scene.quadtree
203 obs = quadtree.leaf_medians
205 if self.misfit_config.optimise_orbital_ramp:
206 stat_level = num.full_like(obs, self.parameter_values['offset'])
208 stat_level += (quadtree.leaf_center_distance[:, 0]
209 * self.parameter_values['ramp_east'])
210 stat_level += (quadtree.leaf_center_distance[:, 1]
211 * self.parameter_values['ramp_north'])
212 statics['displacement.los'] += stat_level
214 stat_syn = statics['displacement.los']
216 res = obs - stat_syn
218 misfit_value = res
219 misfit_norm = obs
221 mf = num.vstack([misfit_value, misfit_norm]).T
222 result = SatelliteMisfitResult(
223 misfits=mf)
225 if self._result_mode == 'full':
226 result.statics_syn = statics
227 result.statics_obs = {'displacement.los': obs}
229 return result
231 def get_combined_weight(self):
232 if self._combined_weight is None:
233 # invcov = self.scene.covariance.weight_matrix
234 # self._combined_weight = invcov * self.manual_weight
235 self._combined_weight = num.full(self.nmisfits, self.manual_weight)
237 return self._combined_weight
239 def prepare_modelling(self, engine, source, targets):
240 return [self]
242 def finalize_modelling(
243 self, engine, source, modelling_targets, modelling_results):
244 return modelling_results[0]
246 def init_bootstrap_residuals(self, nbootstraps, rstate=None, nthreads=0):
247 logger.info(
248 'Scene "%s", initializing bootstrapping residuals from noise '
249 'pertubations...' % self.scene_id)
251 if rstate is None:
252 rstate = num.random.RandomState()
254 scene = self.scene
255 qt = scene.quadtree
256 cov = scene.covariance
257 bootstraps = num.zeros((nbootstraps, qt.nleaves))
259 try:
260 # TODO:mi Signal handler is not given back to the main task!
261 # This is a python3.7 bug
262 logger.warning('Using multi-threading for SatelliteTargets. '
263 'Python 3.7 is buggy and needs to be killed hard:'
264 ' `killall grond`')
265 from concurrent.futures import ThreadPoolExecutor
266 nthreads = os.cpu_count() if not nthreads else nthreads
268 with ThreadPoolExecutor(max_workers=nthreads) as executor:
269 res = executor.map(
270 cov.getQuadtreeNoise,
271 [rstate for _ in range(nbootstraps)])
273 for ibs, bs in enumerate(res):
274 bootstraps[ibs, :] = bs
275 except ImportError:
276 for ibs in range(nbootstraps):
277 if not (ibs+1) % 5:
278 logger.info('Calculating noise realisation %d/%d.'
279 % (ibs+1, nbootstraps))
280 bootstraps[ibs, :] = cov.getQuadtreeNoise(rstate=rstate)
282 self.set_bootstrap_residuals(bootstraps)
284 @classmethod
285 def get_plot_classes(cls):
286 from . import plot
287 plots = super(SatelliteMisfitTarget, cls).get_plot_classes()
288 plots.extend(plot.get_plot_classes())
289 return plots
292__all__ = '''
293 SatelliteTargetGroup
294 SatelliteMisfitConfig
295 SatelliteMisfitTarget
296 SatelliteMisfitResult
297'''.split()