Compare commits

...

47 Commits

Author SHA1 Message Date
asteinbe f38ea955ce simple depth overlap check; pep8 2 years ago
asteinbe af34d1e08d target plot workaround 2 years ago
asteinbe bacd6477c6 waveform targets plotting try workaround 2 years ago
asteinbe 6e1227d8fc minor changes 2 years ago
asteinbe 10418f7461 Update Multisource to new grond, p2 2 years ago
asteinbe 8fb383c565 Merge branch 'master' of https://github.com/pyrocko/grond into multisource 2 years ago
asteinbe 1884b5950e update p1 2 years ago
Andreas Steinberg bd5340b31d
Merge pull request #22 from janrid/multisource 2 years ago
Andreas Steinberg 407d1dbabc
Merge branch 'multisource' into multisource 2 years ago
Ridderbusch 685ab01ac4 fix label 2 years ago
Jan Ridderbusch e59e6aac99 Changed time discretisation of CombiSource problem to optimise time separated sources. 2 years ago
Jan Ridderbusch ac41f9c6b4 Fix satellite closeup plot for multisources. 2 years ago
Jan Ridderbusch 56563ad12e In src/report/app/js/app.js add !grond.MultiRectangularProblem and !grond.DoubleDCProblem as possible yaml structures. 2 years ago
asteinbe 2f968f847d combisource with induvidal times 2 years ago
Jan Ridderbusch ad4111b0a0 In src/problems/multirectangular/problem.py change parameter labels in class MultirectangularProblem. 2 years ago
Jan Ridderbusch e7bd79a470 In src/targets/waveform/plot.py update problem.get_source() for multisource usage. 2 years ago
Jan Ridderbusch 755028761b In src/problems/multirectangular/problem.py change get_source() in class MultirectangularProblem. 2 years ago
miili 6025be0d08 update changelog, bumping patch 2 years ago
miili 05ffcdbbb5 satellite.plot: geographical aspect ratio for LatLon plots 2 years ago
Mi! f5980c6365
Merge pull request #19 from pyrocko/docfix 2 years ago
fabi 8a9d38a5cc docs: corrected files in command for generating reports of examples 2 years ago
fabi 2a808c2185 docs: typo fix 2 years ago
Mi! 5e4e771ebd
Bugfix `GNSSDistributionPlot` 2 years ago
Mi! 2b6b5dc3c8
Merge pull request #18 from pyrocko/issue-template 2 years ago
Mi! 8fa6ef2d48
Update issue templates 2 years ago
Jan Ridderbusch 3e40f3a511 Add/Set nsources to 2. 3 years ago
Jan Ridderbusch 64952ef6b0 Merge branch 'multisource' of https://github.com/janrid/grond into multisource 3 years ago
asteinbe e4d4ad016e nsources None default 3 years ago
braunfuss 097f8fb5fc argv for nsources 3 years ago
braunfuss ceff483ed6 waveform params for multi 3 years ago
asteinbe e3ecde1742 multirectangular source, broken config link of nsources 3 years ago
Jan Ridderbusch 9e8b002a38 Merge branch 'master' into ahar 3 years ago
Jan Ridderbusch 3331bfb3f3 Merge branch 'master' of https://github.com/janrid/grond with pyrocko/grond-master 3 years ago
Jan Ridderbusch 30b647b4f0 Merge branch 'ahar' of https://github.com/janrid/grond into ahar 3 years ago
Jan Ridderbusch 94f3264764 in src/dataset.py if multiple response files per station available now the more recent one will be used. changed mehthod: get_response 3 years ago
Jan Ridderbusch f2c84e1621 Merge branch 'master' of https://github.com/janrid/grond with pyrocko-grond-master 3 years ago
Jan Ridderbusch 72197da305 Merge branch 'master' of janrid/grond into ahar 3 years ago
Jan Ridderbusch 2bce13e403 in src/dataset.py if multiple response files per station available now the more recent one will be used. changed mehthod: get_response 3 years ago
Jan Ridderbusch 993c687155 Merge branch 'master' of https://github.com/janrid/grond 3 years ago
Jan Ridderbusch 00d5114014 in src/dataset.py if multiple response files per station available now the more recent one will be used. changed mehthod: get_response 3 years ago
Jan Ridderbusch 14c7c911fa changed histogram plot default to histogram 3 years ago
Sebastian Heimann 487c0ec919 improve whitelisting implementation 3 years ago
asteinbe 7e3fe6272c nsl station code fix 3 years ago
asteinbe f7bdb851c4 nsources None default 3 years ago
braunfuss 535fbe69ad argv for nsources 3 years ago
braunfuss e04e943d7f waveform params for multi 3 years ago
asteinbe db564822b7 multirectangular source, broken config link of nsources 3 years ago
  1. 20
      .github/ISSUE_TEMPLATE/bug_report.md
  2. 20
      .github/ISSUE_TEMPLATE/feature_request.md
  3. 11
      CHANGELOG.md
  4. 2
      docs/source/examples/satellite_insar/index.rst
  5. 4
      docs/source/examples/waveform_regional/index.rst
  6. 2
      docs/source/examples/waveform_wphase/index.rst
  7. 3
      setup.py
  8. 15
      src/core.py
  9. 20
      src/dataset.py
  10. 266
      src/optimisers/highscore/optimiser.py
  11. 1
      src/problems/__init__.py
  12. 125
      src/problems/base.py
  13. 0
      src/problems/multirectangular/__init__.py
  14. 138
      src/problems/multirectangular/problem.py
  15. 2
      src/report/app/js/app.js
  16. 26
      src/targets/gnss_campaign/plot.py
  17. 176
      src/targets/satellite/plot.py
  18. 44
      src/targets/waveform/plot.py

20
.github/ISSUE_TEMPLATE/bug_report.md

@ -0,0 +1,20 @@
---
name: Bug report
about: Create a report to help us improve
title: ''
labels: bug
assignees: ''
---
**Describe the bug**
A clear and concise description of what the bug is.
**To Reproduce**
Steps to reproduce the behavior:
**Expected behavior**
A clear and concise description of what you expected to happen.
**Additional context**
Add any other context about the problem here.

20
.github/ISSUE_TEMPLATE/feature_request.md

@ -0,0 +1,20 @@
---
name: Feature request
about: Suggest an idea for this project
title: ''
labels: feature request
assignees: ''
---
**Is your feature request related to a problem? Please describe.**
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]
**Describe the solution you'd like**
A clear and concise description of what you want to happen.
**Describe alternatives you've considered**
A clear and concise description of any alternative solutions or features you've considered.
**Additional context**
Add any other context or screenshots about the feature request here.

11
CHANGELOG.md

@ -4,6 +4,17 @@ All notable changes to Grond will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
## [1.0.2]
### Fixed
- Satellite plot: Setting geographical aspect ratio for LatLon data
## [1.0.1]
### Fixed
- GNSS Plotting function
## [1.0.0]
### Added

2
docs/source/examples/satellite_insar/index.rst

@ -170,7 +170,7 @@ Once the optimisation is finished we can generate and open the final report with
.. code-block :: sh
grond report -so rundir/rect_source.grun
grond report -so runs/rect_2009LaAquila.grun
Example report

4
docs/source/examples/waveform_regional/index.rst

@ -10,7 +10,7 @@ To repeat this exercise on your machine, you should first `install Pyrocko <http
.. code-block :: sh
grond init example_regional-cmt grond-playground-regional/
grond init example_regional_cmt grond-playground-regional/
The project folder
@ -180,7 +180,7 @@ Once the optimisation is finished we can generate and open the final report with
.. code-block :: sh
grond report -so config/regional_cmt.gronf
grond report -so runs/cmt_gfz2018pmjk.grun
Example report

2
docs/source/examples/waveform_wphase/index.rst

@ -178,7 +178,7 @@ Once the optimisation is finished we can generate and open the final report with
.. code-block :: sh
grond report -so config/wphase_cmt.gronf
grond report -so runs/wphase_cmt_gfz2015sfdd.grun
Example report
~~~~~~~~~~~~~~

3
setup.py

@ -8,7 +8,7 @@ from setuptools.command.install import install
from setuptools.command.build_py import build_py
version = '1.0.0'
version = '1.0.2'
class CustomInstallCommand(install):
@ -111,6 +111,7 @@ setup(
'grond.problems.cmt',
'grond.problems.double_dc',
'grond.problems.rectangular',
'grond.problems.multirectangular',
'grond.optimisers',
'grond.optimisers.highscore',
'grond.analysers',

15
src/core.py

@ -275,6 +275,7 @@ def check(
results_list = []
sources = []
nsources = 2 # for testing
if n_random_synthetics == 0:
x = problem.get_reference_model()
sources.append(problem.base_source)
@ -284,7 +285,8 @@ def check(
else:
for i in range(n_random_synthetics):
x = problem.get_random_model()
sources.append(problem.get_source(x))
for j in range(nsources): #test
sources.append(problem.get_source(x,j))
results = problem.evaluate(x)
results_list.append(results)
@ -638,21 +640,26 @@ def export(what, rundirs, type=None, pnames=None, filename=None):
print('#', ' '.join(['%16s' % x for x in pnames]), file=out)
def dump(x, gm, indices):
nsources = 2 # fix for testing
sources = []
events = []
if type == 'vector':
print(' ', ' '.join(
'%16.7g' % problem.extract(x, i) for i in indices),
'%16.7g' % gm, file=out)
elif type == 'source':
source = problem.get_source(x)
for i in range(nsources):
source = problem.get_source(x,i)
guts.dump(source, stream=out)
elif type == 'event':
ev = problem.get_source(x).pyrocko_event()
for i in range(nsources):
ev = problem.get_source(x,i).pyrocko_event()
model.dump_events([ev], stream=out)
elif type == 'event-yaml':
ev = problem.get_source(x).pyrocko_event()
ev = problem.get_source(x,0).pyrocko_event()
guts.dump_all([ev], stream=out)
else:

20
src/dataset.py

@ -488,7 +488,21 @@ class Dataset(object):
elif len(candidates) == 0:
raise NotFound('No response found:', (net, sta, loc, cha))
else:
raise NotFound('Multiple responses found:', (net, sta, loc, cha))
# if multiple responses are found take the younger one
t1 = []
index = []
indexSX = -1
for sx in self.responses_stationxml:
indexSX += 1
for network in sx.network_list:
if network.code == net:
# step through all the stations per network
for station in network.station_list:
if station.code == sta:
if t1 == [] or t1 < station.start_date:
t1 = station.start_date
index = indexSX
return candidates[index]
def get_waveform_raw(
self, obj,
@ -710,7 +724,6 @@ class Dataset(object):
for matrix, in_channels, out_channels in projections:
deps = trace.project_dependencies(
matrix, in_channels, out_channels)
try:
trs_restituted_group = []
trs_raw_group = []
@ -741,7 +754,8 @@ class Dataset(object):
trs_raw.extend(trs_raw_group)
except NotFound as e:
exceptions.append((in_channels, out_channels, e))
print((in_channels, out_channels, e))
exceptions.append((in_channels, out_channels, e))
if not trs_projected:
err = []

266
src/optimisers/highscore/optimiser.py

@ -6,9 +6,9 @@ import logging
import time
import numpy as num
from collections import OrderedDict
from pyrocko.guts import StringChoice, Int, Float, Object, List
from pyrocko.guts_array import Array
from shapely.geometry import Polygon, LineString
from grond.meta import GrondError, Forbidden, has_get_plot_classes
from grond.problems.base import ModelHistory
@ -18,7 +18,7 @@ from grond.optimisers.base import Optimiser, OptimiserConfig, BadProblem, \
guts_prefix = 'grond'
logger = logging.getLogger('grond.optimisers.highscore.optimiser')
d2r = math.pi / 180.
def nextpow2(i):
return 2**int(math.ceil(math.log(i)/math.log(2.)))
@ -142,9 +142,7 @@ class SamplerPhase(Object):
class InjectionSamplerPhase(SamplerPhase):
xs_inject = Array.T(
dtype=num.float, shape=(None, None),
help='Array with the reference model.')
xs_inject = Array.T(dtype=num.float, shape=(None, None))
def get_raw_sample(self, problem, iiter, chains):
return Sample(model=self.xs_inject[iiter, :])
@ -154,28 +152,55 @@ class UniformSamplerPhase(SamplerPhase):
def get_raw_sample(self, problem, iiter, chains):
xbounds = problem.get_parameter_bounds()
intersect = True
while intersect is True:
pars = problem.random_uniform(xbounds, self.get_rstate())
nsources = 2
sources = []
if nsources is not None:
for i in range(nsources):
source = problem.get_source(pars, i)
sources.append(source)
depths_max = []
depths_min = []
polygons = []
for src in sources:
depths_max.append(num.max(src.outline()[:, 2]))
depths_min.append(num.min(src.outline()[:, 2]))
src_outline = sources[0].outline('xy')
polygons.append(Polygon(src_outline))
for k in range(len(polygons)):
for j in range(len(polygons)):
p1 = polygons[k]
p2 = polygons[j]
if not p1.intersects(p2) or p1 == p2:
intersect = False
else:
line_1 = [(1.0, depths_min[k]), (1.0, depths_max[k])]
line_2 = [(1.0, depths_min[j]), (1.0, depths_max[j])]
line1 = LineString(line_1)
line2 = LineString(line_2)
intersection_depth = line1.intersection(line2)
if intersection_depth is False:
intersect = False
else:
intersect = True
print('intersection, uniform phase redraw')
if any(sources.count(x) > 1 for x in sources):
intersect = True
return Sample(model=problem.random_uniform(xbounds, self.get_rstate()))
class DirectedSamplerPhase(SamplerPhase):
scatter_scale = Float.T(
optional=True,
help='Scales search radius around the current `highscore` models')
scatter_scale_begin = Float.T(
optional=True,
help='Scaling factor at beginning of the phase.')
scatter_scale_end = Float.T(
optional=True,
help='Scaling factor at the end of the directed phase.')
scatter_scale = Float.T(optional=True)
scatter_scale_begin = Float.T(optional=True)
scatter_scale_end = Float.T(optional=True)
starting_point = SamplerStartingPointChoice.T(
default='excentricity_compensated',
help='Tunes to the center value of the sampler distribution.'
'May increase the likelihood to draw a highscore member model'
' off-center to the mean value')
default='excentricity_compensated')
sampler_distribution = SamplerDistributionChoice.T(
default='normal',
help='Distribution new models are drawn from.')
default='normal')
standard_deviation_estimator = StandardDeviationEstimatorChoice.T(
default='median_density_single_chain')
@ -208,36 +233,63 @@ class DirectedSamplerPhase(SamplerPhase):
ilink_choice = None
ichain_choice = num.argmin(chains.accept_sum)
intersect = True
while intersect is True:
if self.starting_point == 'excentricity_compensated':
models = chains.models(ichain_choice)
ilink_choice = excentricity_compensated_choice(
models,
chains.standard_deviation_models(
ichain_choice, self.standard_deviation_estimator),
2., rstate)
if self.starting_point == 'excentricity_compensated':
models = chains.models(ichain_choice)
ilink_choice = excentricity_compensated_choice(
models,
chains.standard_deviation_models(
ichain_choice, self.standard_deviation_estimator),
2., rstate)
xchoice = chains.model(ichain_choice, ilink_choice)
xchoice = chains.model(ichain_choice, ilink_choice)
elif self.starting_point == 'random':
ilink_choice = rstate.randint(0, chains.nlinks)
xchoice = chains.model(ichain_choice, ilink_choice)
elif self.starting_point == 'random':
ilink_choice = rstate.randint(0, chains.nlinks)
xchoice = chains.model(ichain_choice, ilink_choice)
elif self.starting_point == 'mean':
xchoice = chains.mean_model(ichain_choice)
else:
assert False, 'invalid starting_point choice: %s' % (
self.starting_point)
elif self.starting_point == 'mean':
xchoice = chains.mean_model(ichain_choice)
ntries_sample = 0
if self.sampler_distribution == 'normal':
x = num.zeros(npar, dtype=num.float)
sx = chains.standard_deviation_models(
ichain_choice, self.standard_deviation_estimator)
for ipar in range(npar):
ntries = 0
else:
assert False, 'invalid starting_point choice: %s' % (
self.starting_point)
ntries_sample = 0
if self.sampler_distribution == 'normal':
x = num.zeros(npar, dtype=num.float)
sx = chains.standard_deviation_models(
ichain_choice, self.standard_deviation_estimator)
for ipar in range(npar):
ntries = 0
while True:
if sx[ipar] > 0.:
v = num.random.normal(
xchoice[ipar],
factor*sx[ipar])
else:
v = xchoice[ipar]
if xbounds[ipar, 0] <= v and \
v <= xbounds[ipar, 1]:
break
if ntries > self.ntries_sample_limit:
raise GrondError(
'failed to produce a suitable '
'candidate sample from normal '
'distribution')
ntries += 1
x[ipar] = v
elif self.sampler_distribution == 'multivariate_normal':
ok_mask_sum = num.zeros(npar, dtype=num.int)
while True:
if sx[ipar] > 0.:
v = rstate.normal(
@ -249,50 +301,86 @@ class DirectedSamplerPhase(SamplerPhase):
if xbounds[ipar, 0] <= v and \
v <= xbounds[ipar, 1]:
break
ok_mask = num.logical_and(
xbounds[:, 0] <= xcandi, xcandi <= xbounds[:, 1])
if ntries > self.ntries_sample_limit:
logger.warning(
'failed to produce a suitable '
'candidate sample from normal '
'distribution for parameter \'%s\''
'- drawing from uniform instead.' %
pnames[ipar])
v = rstate.uniform(xbounds[ipar, 0],
xbounds[ipar, 1])
break
if ntries > self.ntries_sample_limit:
logger.warning(
'failed to produce a suitable '
'candidate sample from normal '
'distribution for parameter \'%s\''
'- drawing from uniform instead.' %
pnames[ipar])
v = rstate.uniform(xbounds[ipar, 0],
xbounds[ipar, 1])
break
ntries += 1
ntries += 1
x[ipar] = v
x[ipar] = v
elif self.sampler_distribution == 'multivariate_normal':
ok_mask_sum = num.zeros(npar, dtype=num.int)
while True:
ntries_sample += 1
xcandi = rstate.multivariate_normal(
xchoice, factor**2 * chains.cov(ichain_choice))
elif self.sampler_distribution == 'multivariate_normal':
ok_mask_sum = num.zeros(npar, dtype=num.int)
while True:
ntries_sample += 1
xcandi = rstate.multivariate_normal(
xchoice, factor**2 * chains.cov(ichain_choice))
ok_mask = num.logical_and(
xbounds[:, 0] <= xcandi, xcandi <= xbounds[:, 1])
ok_mask = num.logical_and(
xbounds[:, 0] <= xcandi, xcandi <= xbounds[:, 1])
if num.all(ok_mask):
break
if num.all(ok_mask):
break
ok_mask_sum += ok_mask
ok_mask_sum += ok_mask
if ntries_sample > self.ntries_sample_limit:
logger.warning(
'failed to produce a suitable candidate '
'sample from multivariate normal '
'distribution, (%s) - drawing from uniform instead' %
', '.join('%s:%i' % xx for xx in
zip(pnames, ok_mask_sum)))
xbounds = problem.get_parameter_bounds()
xcandi = problem.random_uniform(xbounds, rstate)
break
if ntries_sample > self.ntries_sample_limit:
logger.warning(
'failed to produce a suitable candidate '
'sample from multivariate normal '
'distribution, (%s) - drawing from uniform instead' %
', '.join('%s:%i' % xx for xx in
zip(pnames, ok_mask_sum)))
xbounds = problem.get_parameter_bounds()
xcandi = problem.random_uniform(xbounds, rstate)
break
x = xcandi
x = xcandi
pars = x
nsources = 2
sources = []
if nsources is not None:
for i in range(nsources):
source = problem.get_source(pars, i)
sources.append(source)
depths_max = []
depths_min = []
polygons = []
for src in sources:
depths_max.append(num.max(src.outline()[:,2]))
depths_min.append(num.min(src.outline()[:,2]))
src_outline = sources[0].outline('xy')
polygons.append(Polygon(src_outline))
for k in range(len(polygons)):
for j in range(len(polygons)):
p1 = polygons[k]
p2 = polygons[j]
if not p1.intersects(p2) or p1 == p2:
intersect = False
else:
line_1 = [(1.0, depths_min[k]), (1.0, depths_max[k])]
line_2 = [(1.0, depths_min[j]), (1.0, depths_max[j])]
line1 = LineString(line_1)
line2 = LineString(line_2)
intersection_depth = line1.intersection(line2)
if intersection_depth is False:
intersect = False
else:
intersect = True
print('intersection, uniform phase redraw')
if any(sources.count(x) > 1 for x in sources):
intersect = True
return Sample(
model=x,
@ -333,9 +421,9 @@ class Chains(object):
self.nchains = nchains
self.nlinks_cap = nlinks_cap
self.chains_m = num.zeros(
(self.nchains, nlinks_cap), dtype=num.float)
(self.nchains, nlinks_cap), num.float)
self.chains_i = num.zeros(
(self.nchains, nlinks_cap), dtype=num.int)
(self.nchains, nlinks_cap), num.int)
self.nlinks = 0
self.nread = 0
@ -738,17 +826,9 @@ class HighScoreOptimiserConfig(OptimiserConfig):
sampler_phases = List.T(
SamplerPhase.T(),
default=[UniformSamplerPhase(niterations=1000),
DirectedSamplerPhase(niterations=5000)],
help='Stages of the sampler: Start with uniform sampling of the model'
' model space and narrow down through directed sampling.')
chain_length_factor = Float.T(
default=8.,
help='Controls the length of each chain: '
'chain_length_factor * nparameters + 1')
nbootstrap = Int.T(
default=100,
help='Number of bootstrap realisations to be tracked simultaneously in'
' the optimisation.')
DirectedSamplerPhase(niterations=5000)])
chain_length_factor = Float.T(default=8.)
nbootstrap = Int.T(default=100)
def get_optimiser(self):
return HighScoreOptimiser(

1
src/problems/__init__.py

@ -5,4 +5,5 @@ Definition of the objective function and source model parameter space.
from .base import * # noqa
from .cmt.problem import * # noqa
from .rectangular.problem import * # noqa
from .multirectangular.problem import * # noqa
from .double_dc.problem import * # noqa

125
src/problems/base.py

@ -25,7 +25,6 @@ from grond import stats
from grond.version import __version__
guts_prefix = 'grond'
logger = logging.getLogger('grond.problems.base')
km = 1e3
@ -33,11 +32,64 @@ as_km = dict(scale_factor=km, scale_unit='km')
g_rstate = num.random.RandomState()
def range_overlap(a_min, a_max, b_min, b_max):
'''Neither range is completely greater than the other
'''
overlapping = True
if (a_min > b_max) or (a_max < b_min):
overlapping = False
return overlapping
def overlap(r1, r2):
'''Overlapping rectangles overlap both horizontally & vertically
'''
return range_overlap(r1.left, r1.right, r2.left, r2.right) and range_overlap(r1.bottom, r1.top, r2.bottom, r2.top)
def nextpow2(i):
return 2**int(math.ceil(math.log(i)/math.log(2.)))
class CombiSource(gf.Source):
'''Composite source model.'''
discretized_source_class = gf.DiscretizedMTSource
subsources = List.T(gf.Source.T())
def __init__(self, subsources=[], **kwargs):
if subsources:
lats = num.array(
[subsource.lat for subsource in subsources], dtype=num.float)
lons = num.array(
[subsource.lon for subsource in subsources], dtype=num.float)
assert num.all(lats == lats[0]) and num.all(lons == lons[0])
lat, lon = lats[0], lons[0]
depth = float(num.mean([p.depth for p in subsources]))
t = float(list([p.time for p in subsources])[0])
kwargs.update(time=t, lat=float(lat), lon=float(lon), depth=depth)
gf.Source.__init__(self, subsources=subsources, **kwargs)
def get_factor(self):
return 1.0
def discretize_basesource(self, store, target=None):
dsources = []
t0 = self.subsources[0].time
t1 = self.subsources[1].time
tdiff = t0-t1
for sf in self.subsources:
ds = sf.discretize_basesource(store, target)
ds.m6s *= sf.get_factor()
dsources.append(ds)
dsources[1].times = dsources[1].times - tdiff
return gf.DiscretizedMTSource.combine(dsources)
class ProblemConfig(Object):
'''
Base class for config section defining the objective function setup.
@ -45,9 +97,11 @@ class ProblemConfig(Object):
Factory for :py:class:`Problem` objects.
'''
name_template = String.T()
nsources = Int.T(default=1)
norm_exponent = Int.T(default=2)
nthreads = Int.T(default=1)
def get_problem(self, event, target_groups, targets):
'''
Instantiate the problem with a given event and targets.
@ -69,17 +123,16 @@ class Problem(Object):
dependants = List.T(Parameter.T())
norm_exponent = Int.T(default=2)
base_source = gf.Source.T(optional=True)
nsources = Int.T(default=1)
targets = List.T(MisfitTarget.T())
target_groups = List.T(TargetGroup.T())
grond_version = String.T(optional=True)
nthreads = Int.T(default=1)
def __init__(self, **kwargs):
Object.__init__(self, **kwargs)
if self.grond_version is None:
self.grond_version = __version__
self._target_weights = None
self._engine = None
self._family_mask = None
@ -89,6 +142,7 @@ class Problem(Object):
self.problem_parameters + self.problem_waveform_parameters
self.check()
logger.name = self.__class__.__name__
@classmethod
def get_plot_classes(cls):
@ -106,9 +160,6 @@ class Problem(Object):
paths.add(grp.path)
logger.debug('TargetGroup check OK.')
def get_engine(self):
return self._engine
def copy(self):
o = copy.copy(self)
o._target_weights = None
@ -120,11 +171,19 @@ class Problem(Object):
target.set_parameter_values(x[nprob:nprob+target.nparameters])
nprob += target.nparameters
def get_parameter_dict(self, model, group=None):
def get_parameter_dict(self, model, group=None, nsources=None):
params = []
for ip, p in enumerate(self.parameters):
if group in p.groups or group is None:
params.append((p.name, model[ip]))
if nsources is not None:
for ip, p in enumerate(self.parameters):
if group in p.groups:
try:
params.append((p.name[:-1], model[ip]))
except:
pass
else:
for ip, p in enumerate(self.parameters):
if group in p.groups:
params.append((p.name, model[ip]))
return ADict(params)
def get_parameter_array(self, d):
@ -251,6 +310,14 @@ class Problem(Object):
def set_engine(self, engine):
self._engine = engine
def get_engine(self):
return self._engine
def get_gf_store(self, target):
if self.get_engine() is None:
raise GrondError('Cannot get GF Store, modelling is not set up!')
return self.get_engine().get_store(target.store_id)
def random_uniform(self, xbounds, rstate):
x = rstate.uniform(0., 1., self.nparameters)
x *= (xbounds[:, 1] - xbounds[:, 0])
@ -467,10 +534,21 @@ class Problem(Object):
return self._family_mask
def evaluate(self, x, mask=None, result_mode='full', targets=None):
source = self.get_source(x)
engine = self.get_engine()
def evaluate(self, x, mask=None, result_mode='full',
targets=None, nsources=None):
patches = []
nsources = 2
outlines = []
if nsources == 2:
for i in range(nsources):
source = self.get_source(x, i)
patches.append(source)
outlines.append(source.outline())
else:
source = self.get_source(x)
engine = self.get_engine()
sources = CombiSource(subsources=patches)
self.set_target_parameter_values(x)
if mask is not None and targets is not None:
@ -483,7 +561,9 @@ class Problem(Object):
modelling_targets = []
t2m_map = {}
for itarget, target in enumerate(targets):
t2m_map[target] = target.prepare_modelling(engine, source, targets)
t2m_map[target] = target.prepare_modelling(engine,
sources,
targets)
if mask is None or mask[itarget]:
modelling_targets.extend(t2m_map[target])
@ -496,8 +576,9 @@ class Problem(Object):
modelling_targets_unique = list(u2m_map.keys())
resp = engine.process(source, modelling_targets_unique,
resp = engine.process(sources, modelling_targets_unique,
nthreads=self.nthreads)
modelling_results_unique = list(resp.results_list[0])
modelling_results = [None] * len(modelling_targets)
@ -514,14 +595,14 @@ class Problem(Object):
nmt_this = len(t2m_map[target])
if mask is None or mask[itarget]:
result = target.finalize_modelling(
engine, source,
engine, sources,
t2m_map[target],
modelling_results[imt:imt+nmt_this])
imt += nmt_this
else:
result = gf.SeismosizerError(
'target was excluded from modelling')
'target was excluded from modelling')
results.append(result)
@ -560,10 +641,10 @@ class Problem(Object):
return results
def get_random_model(self):
def get_random_model(self, ntries_limit=100):
xbounds = self.get_parameter_bounds()
while True:
for _ in range(ntries_limit):
x = self.random_uniform(xbounds, rstate=g_rstate)
try:
return self.preconstrain(x)
@ -571,6 +652,10 @@ class Problem(Object):
except Forbidden:
pass
raise GrondError(
'Could not find any suitable candidate sample within %i tries' % (
ntries_limit))
class ProblemInfoNotAvailable(GrondError):
pass

0
src/problems/multirectangular/__init__.py

138
src/problems/multirectangular/problem.py

@ -0,0 +1,138 @@
import numpy as num
import logging
import sys
from pyrocko import gf, util
from pyrocko.guts import String, Float, Dict, Int
from grond.meta import expand_template, Parameter
from ..base import Problem, ProblemConfig
guts_prefix = 'grond'
logger = logging.getLogger('grond.problems.multirectangular.problem')
km = 1e3
as_km = dict(scale_factor=km, scale_unit='km')
class MultiRectangularProblemConfig(ProblemConfig):
ranges = Dict.T(String.T(), gf.Range.T())
decimation_factor = Int.T(default=1)
distance_min = Float.T(default=0.)
nsources = Int.T(default=1)
def get_problem(self, event, target_groups, targets):
base_source = gf.RectangularSource.from_pyrocko_event(
event,
anchor='top',
decimation_factor=self.decimation_factor)
subs = dict(
event_name=event.name,
event_time=util.time_to_str(event.time))
problem = MultiRectangularProblem(
name=expand_template(self.name_template, subs),
base_source=base_source,
distance_min=self.distance_min,
target_groups=target_groups,
targets=targets,
ranges=self.ranges,
norm_exponent=self.norm_exponent)
return problem
class MultiRectangularProblem(Problem):
nsources = 2
for i in range(0, 100):
if "--nsources="+str(i) in sys.argv:
nsources = int(i)
if nsources is None:
print('input --nsources= to go command missing')
problem_parameters = []
problem_waveform_parameters = []
for i in range(0, nsources):
problem_parameters.append(Parameter('north_shift%s' % i,
'm',
label='Northing %s' %i,
**as_km))
problem_parameters.append(Parameter('east_shift%s' % i,
'm',
label='Easting %s' %i,
**as_km))
problem_parameters.append(Parameter('depth%s' % i, 'm',
label='Depth %s' %i,
**as_km))
problem_parameters.append(Parameter('length%s' % i,
'm',
label='Length %s' %i,
**as_km))
problem_parameters.append(Parameter('width%s' % i, 'm',
label='Width %s' %i,
**as_km))
problem_parameters.append(Parameter('dip%s' % i, 'deg',
label='Dip %s' %i))
problem_parameters.append(Parameter('strike%s' % i,
'deg',
label='Strike %s' %i))
problem_parameters.append(Parameter('rake%s' % i, 'deg',
label='Rake %s' %i))
problem_parameters.append(Parameter('slip%s' % i, 'm',
label='Slip %s' %i))
problem_parameters.append(Parameter('nucleation_x%s' % i,
'offset',
label='Nucleation X %s' %i))
problem_parameters.append(Parameter('nucleation_y%s' % i,
'offset',
label='Nucleation Y %s' %i))
problem_parameters.append(Parameter('time%s' % i, 's',
label='Time %s' %i))
dependants = []
distance_min = Float.T(default=0.0)
def pack(self, source):
arr = self.get_parameter_array(source)
for ip, p in enumerate(self.parameters):
if p.name == 'time':
arr[ip] -= self.base_source.time
return arr
def get_source(self, x, i):
d = self.get_parameter_dict(x[0+12*i:12+i*12], nsources=True)
p = {}
for k in self.base_source.keys():
if k in d:
p[k] = float(
self.ranges[k+str(i)].make_relative(self.base_source[k],
d[k]))
source = self.base_source.clone(**p)
return source
def random_uniform(self, xbounds, rstate):
x = num.zeros(self.nparameters)
for i in range(self.nparameters):
x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
return x
def preconstrain(self, x):
# source = self.get_source(x)
# if any(self.distance_min > source.distance_to(t)
# for t in self.targets):
# raise Forbidden()
return x
@classmethod
def get_plot_classes(cls):
plots = super(MultiRectangularProblem, cls).get_plot_classes()
return plots
__all__ = '''
MultiRectangularProblem
MultiRectangularProblemConfig
'''.split()

2
src/report/app/js/app.js

@ -115,6 +115,8 @@ var yaml_type_map = [
['!grond.FeatureMeasure', Dummy],
['!grond.CMTProblem', Dummy],
['!grond.RectangularProblem', Dummy],
['!grond.DoubleDCProblem', Dummy],
['!grond.MultiRectangularProblem', Dummy],
['!pf.MTSource', Dummy],
['!pf.RectangularSource', Dummy],
['!pf.HalfSinusoidSTF', Dummy],

26
src/targets/gnss_campaign/plot.py

@ -7,7 +7,7 @@ from pyrocko import orthodrome as od
from grond.plot.config import PlotConfig
from grond.plot.collection import PlotItem
from grond.problems import CMTProblem, RectangularProblem
from grond.problems import CMTProblem, RectangularProblem, MultiRectangularProblem
from ..plot import StationDistributionPlot
@ -79,8 +79,15 @@ the upper fault edge.
gms = gms[isort]
models = history.models[isort, :]
xbest = models[0, :]
source = problem.get_source(xbest)
nsources = 2
if nsources is not None:
sources = []
for i in range(nsources):
source = problem.get_source(xbest, i)
sources.append(source)
else:
source = problem.get_source(xbest)
results = problem.evaluate(
xbest, result_mode='full', targets=gnss_targets)
@ -189,6 +196,17 @@ displacements derived from best rupture model (red).
t=60,
*m.jxyr)
elif isinstance(problem, MultiRectangularProblem):
if nsources is not None:
for subsource in sources:
m.gmt.psxy(
in_rows=subsource.outline(cs='lonlat'),
L='+p2p,black',
W='1p,black',
G='black',
t=60,
*m.jxyr)
return (item, m)
ifig = 0
@ -244,8 +262,6 @@ components).
item = PlotItem(name='station_distribution-%s' % target.path)
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, ws, labels)
fig.suptitle('GNSS Station Distribution and Weight (%s)'
% target.path, fontsize=self.font_size_title)
legend.set_title('Weight')
yield (item, fig)

176
src/targets/satellite/plot.py

@ -1,3 +1,4 @@
import logging
import numpy as num
from matplotlib import cm, gridspec
@ -8,7 +9,10 @@ from matplotlib import pyplot as plt
from matplotlib import patches
from pyrocko.guts import Tuple, Float, String, Int, Bool
km = 1000.
logger = logging.getLogger('grond.targets.satellite.plot')
km = 1e3
d2r = num.pi/180.
guts_prefix = 'grond'
@ -54,20 +58,21 @@ class SatelliteTargetDisplacement(PlotConfig):
title=u'InSAR Displacements',
section='fits',
feather_icon='navigation',
description=u'''
Maps showing subsampled surface displacements as observed, modelled and the
residual (observed minus modelled).
The displacement values predicted by the orbit-ambiguity ramps are added to the
modelled displacements (middle panels). The color shows the LOS displacement
values associated with, and the extent of, every quadtree box. The light grey
dots show the focal point of pixels combined in the quadtree box. This point
corresponds to the position of the modelled data point.
The large dark grey dot shows the reference source position. The grey filled
box shows the surface projection of the modelled source, with the thick-lined
edge marking the upper fault edge. Complete data extent is shown.
''')
description=u' Maps showing subsampled surface displacements as'
u' observed, modelled and the residual (observed minus'
u' modelled).\n The displacement values predicted by'
u' the orbit-ambiguity ramps are added to the modelled'
u' displacements (middle panels). The color shows the'
u' LOS displacement values associated with, and the'
u' extent of, every quadtree box. The light grey dots'
u' show the focal point of pixels combined in the'
u' quadtree box. This point corresponds to the'
u' position of the modelled data point.\n The large dark'
u' grey dot shows the reference source position. The'
u' grey filled box shows the surface projection of the'
u' modelled source, with the thick-lined edge marking'
u' the upper fault edge. '
u' Complete data extent is shown.')
def draw_static_fits(self, ds, history, optimiser, closeup=False):
from pyrocko.orthodrome import latlon_to_ne_numpy
@ -82,8 +87,16 @@ edge marking the upper fault edge. Complete data extent is shown.
gms = gms[isort]
models = history.models[isort, :]
xbest = models[0, :]
source = problem.get_source(xbest)
# nsources = problem.nsources #help
nsources = 2
if nsources is not None:
sources = []
for i in range(nsources):
source_i = problem.get_source(xbest, i)
sources.append(source_i)
source = sources[0]
else:
source = problem.get_source(xbest)
results = problem.evaluate(xbest, targets=sat_targets)
def initAxes(ax, scene, title, last_axes=False):
@ -97,7 +110,8 @@ edge marking the upper fault edge. Complete data extent is shown.
if not self.relative_coordinates:
import utm
utm_E, utm_N, utm_zone, utm_zone_letter =\
utm.from_latlon(source.lat, source.lon)
utm.from_latlon(source.effective_lat,
source.effective_lon)
scale_x['offset'] = utm_E
scale_y['offset'] = utm_N
@ -107,33 +121,65 @@ edge marking the upper fault edge. Complete data extent is shown.
va='bottom', ha='right',
fontsize=8, alpha=.7,
transform=ax.transAxes)
ax.set_aspect('equal')
elif scene.frame.isDegree():
ax.set_xlabel('Lon [°]')
scale_x = {'scale': 1.}
scale_y = {'scale': 1.}
if not self.relative_coordinates:
scale_x['offset'] = source.lon
scale_y['offset'] = source.lat
scale_x['offset'] = source.effective_lon
scale_y['offset'] = source.effective_lat
ax.set_aspect(1./num.cos(source.effective_lon*d2r))
scale_axes(ax.get_xaxis(), **scale_x)
scale_axes(ax.get_yaxis(), **scale_y)
ax.set_aspect('equal')
def drawSource(ax, scene):
if scene.frame.isMeter():
fn, fe = source.outline(cs='xy').T
if nsources is not None:
fns = []
fes = []
for subsource in sources:
fn_sub, fe_sub = subsource.outline(cs='xy').T
fes.append(fe_sub)
fns.append(fn_sub)
else:
fn, fe = source.outline(cs='xy').T
elif scene.frame.isDegree():
fn, fe = source.outline(cs='latlon').T
fn -= source.lat
fe -= source.lon
if nsources is not None:
fns = []
fes = []
for subsource in sources:
fn_sub, fe_sub = subsource.outline(cs='latlon').T
fn_sub -= subsource.lat
fe_sub -= subsource.lon
fes.append(fe_sub)
fns.append(fn_sub)
else:
fn, fe = source.outline(cs='latlon').T
fn -= source.lat
fe -= source.lon
fn -= source.effective_lat
fe -= source.effective_lon
# source is centered
ax.scatter(0., 0., color='black', s=3, alpha=.5, marker='o')
ax.fill(fe, fn,
edgecolor=(0., 0., 0.),
facecolor=(.5, .5, .5), alpha=0.7)
ax.plot(fe[0:2], fn[0:2], 'k', linewidth=1.3)
if nsources is not None:
for fe, fn in zip(fes, fns):
ax.fill(fe, fn,
edgecolor=(0., 0., 0.),
facecolor=(.5, .5, .5), alpha=0.5)
ax.plot(fe[0:2], fn[0:2], 'k', linewidth=1.3)
else:
ax.fill(fe, fn,
edgecolor=(0., 0., 0.),
facecolor=(.1, .2, .4), alpha=0.5)
ax.plot(fe[0:2], fn[0:2], 'k', linewidth=1.3)
def mapDisplacementGrid(displacements, scene):
arr = num.full_like(scene.displacement, fill_value=num.nan)
@ -197,10 +243,10 @@ edge marking the upper fault edge. Complete data extent is shown.
if target.scene.frame.isMeter():
off_n, off_e = map(float, latlon_to_ne_numpy(
target.scene.frame.llLat, target.scene.frame.llLon,
source.lat, source.lon))
source.effective_lat, source.effective_lon))
if target.scene.frame.isDegree():
off_n = source.lat - target.scene.frame.llLat
off_e = source.lon - target.scene.frame.llLon
off_n = source.effective_lat - target.scene.frame.llLat
off_e = source.effective_lon - target.scene.frame.llLon
turE, turN, tllE, tllN = zip(
*[(l.gridE.max()-off_e,
@ -244,10 +290,10 @@ data and (right) the model residual.
if scene.frame.isMeter():
offset_n, offset_e = map(float, latlon_to_ne_numpy(
scene.frame.llLat, scene.frame.llLon,
source.lat, source.lon))
source.effective_lat, source.effective_lon))
elif scene.frame.isDegree():
offset_n = source.lat - scene.frame.llLat
offset_e = source.lon - scene.frame.llLon
offset_n = source.effective_lat - scene.frame.llLat
offset_e = source.effective_lon - scene.frame.llLon
im_extent = (scene.frame.E.min() - offset_e,
scene.frame.E.max() - offset_e,
@ -312,17 +358,57 @@ data and (right) the model residual.
if closeup:
if scene.frame.isMeter():
fn, fe = source.outline(cs='xy').T
elif scene.frame.isDegree():
fn, fe = source.outline(cs='latlon').T
fn -= source.lat
fe -= source.lon
off_n = (fn[0] + fn[1]) / 2
off_e = (fe[0] + fe[1]) / 2
if nsources is not None:
fns = []
fes = []
for subsource in sources:
fn_sub, fe_sub = subsource.outline(cs='xy').T
fes.append(fe_sub)
fns.append(fn_sub)
fnv = list(fns[0])
fnv.extend(list(fns[1]))
fnv = num.array(fnv)
fev = list(fes[0])
fev.extend(list(fes[1]))
fev = num.array(fev)
off_n = (fns[0][0] + fns[1][1]) / 2
off_e = (fes[0][0] + fes[1][1]) / 2
else:
fn, fe = source.outline(cs='xy').T
off_n = (fn[0] + fn[1]) / 2
off_e = (fe[0] + fe[1]) / 2
fnv = fn
fev = fe
fault_size = 2*num.sqrt(max(abs(fn-off_n))**2
+ max(abs(fe-off_e))**2)
elif scene.frame.isDegree():
if nsources is not None:
fns = []
fes = []
for subsource in sources:
fn_sub, fe_sub = subsource.outline(cs='latlon').T
fn_sub -= subsource.lat
fe_sub -= subsource.lon
fes.append(fe_sub)
fns.append(fn_sub)
fnv = list(fns[0])
fnv.extend(list(fns[1]))
fnv = num.array(fnv)
fev = list(fes[0])
fev.extend