1import logging
2import numpy as num
3from pyrocko.model import gnss
5from pyrocko.plot import automap, mpl_init
6from pyrocko import orthodrome as od
8from grond.plot.config import PlotConfig
9from grond.plot.collection import PlotItem
10from grond.problems import CMTProblem, RectangularProblem, \
11 VLVDProblem
13from ..plot import StationDistributionPlot
15import copy
16from pyrocko.guts import Tuple, Float, Bool
18guts_prefix = 'grond'
19km = 1e3
21logger = logging.getLogger('grond.targets.gnss_campaign.plot')
24class GNSSTargetMisfitPlot(PlotConfig):
25 ''' Maps showing horizontal surface displacements
26 of a GNSS campaign and model '''
28 name = 'gnss'
30 size_cm = Tuple.T(
31 2, Float.T(),
32 default=(30., 30.),
33 help='width and length of the figure in cm')
34 show_topo = Bool.T(
35 default=False,
36 help='show topography')
37 show_grid = Bool.T(
38 default=True,
39 help='show the lat/lon grid')
40 show_rivers = Bool.T(
41 default=True,
42 help='show rivers on the map')
43 radius = Float.T(
44 optional=True,
45 help='radius of the map around campaign center lat/lon')
47 def make(self, environ):
48 cm = environ.get_plot_collection_manager()
49 history = environ.get_history(subset='harvest')
50 optimiser = environ.get_optimiser()
51 ds = environ.get_dataset()
53 environ.setup_modelling()
55 cm.create_group_automap(
56 self,
57 self.draw_gnss_fits(ds, history, optimiser),
58 title=u'GNSS Displacements',
59 section='fits',
60 feather_icon='map',
61 description=u'''
62Maps showing station positions and statiom names of the GNSS targets.
64Arrows the observed surface displacements (black arrows) and synthetic
65displacements (red arrows). The top plot shows the horizontal displacements and
66the bottom plot the vertical displacements. The grey filled box shows the
67surface projection of the modelled source, with the thick-lined edge marking
68the upper fault edge.
69''')
71 def draw_gnss_fits(self, ds, history, optimiser, vertical=False):
72 problem = history.problem
74 gnss_targets = problem.gnss_targets
75 for target in gnss_targets:
76 target.set_dataset(ds)
78 xbest = history.get_best_model()
79 source = history.get_best_source()
81 results = problem.evaluate(
82 xbest, result_mode='full', targets=gnss_targets)
84 def plot_gnss(gnss_target, result, ifig, vertical=False):
85 campaign = gnss_target.campaign
86 item = PlotItem(
87 name='fig_%i' % ifig,
88 attributes={
89 'targets': gnss_target.path
90 },
91 title=u'Static GNSS Surface Displacements - Campaign %s'
92 % campaign.name,
93 description=u'''
94Static surface displacement from GNSS campaign %s (black vectors) and
95displacements derived from best model (red).
96''' % campaign.name)
98 event = source.pyrocko_event()
99 locations = campaign.stations + [event]
101 lat, lon = od.geographic_midpoint_locations(locations)
103 if self.radius is None:
104 coords = num.array([loc.effective_latlon for loc in locations])
105 radius = od.distance_accurate50m_numpy(
106 lat[num.newaxis], lon[num.newaxis],
107 coords[:, 0].max(), coords[:, 1]).max()
108 radius *= 1.1
110 if radius < 30.*km:
111 logger.warning(
112 'Radius of GNSS campaign %s too small, defaulting'
113 ' to 30 km' % campaign.name)
114 radius = 30*km
116 model_camp = gnss.GNSSCampaign(
117 stations=copy.deepcopy(campaign.stations),
118 name='grond model')
119 for ista, sta in enumerate(model_camp.stations):
120 sta.north.shift = result.statics_syn['displacement.n'][ista]
121 sta.north.sigma = 0.
123 sta.east.shift = result.statics_syn['displacement.e'][ista]
124 sta.east.sigma = 0.
126 if sta.up:
127 sta.up.shift = -result.statics_syn['displacement.d'][ista]
128 sta.up.sigma = 0.
130 m = automap.Map(
131 width=self.size_cm[0],
132 height=self.size_cm[1],
133 lat=lat,
134 lon=lon,
135 radius=radius,
136 show_topo=self.show_topo,
137 show_grid=self.show_grid,
138 show_rivers=self.show_rivers,
139 color_wet=(216, 242, 254),
140 color_dry=(238, 236, 230))
142 all_stations = campaign.stations + model_camp.stations
143 offset_scale = num.zeros(len(all_stations))
145 for ista, sta in enumerate(all_stations):
146 for comp in sta.components.values():
147 offset_scale[ista] += comp.shift
148 offset_scale = num.sqrt(offset_scale**2).max()
150 m.add_gnss_campaign(
151 campaign,
152 psxy_style={
153 'G': 'black',
154 'W': '0.8p,black',
155 },
156 offset_scale=offset_scale,
157 vertical=vertical)
159 m.add_gnss_campaign(
160 model_camp,
161 psxy_style={
162 'G': 'red',
163 'W': '0.8p,red',
164 't': 30,
165 },
166 offset_scale=offset_scale,
167 vertical=vertical,
168 labels=False)
170 if isinstance(problem, CMTProblem):
171 from pyrocko import moment_tensor
172 from pyrocko.plot import gmtpy
174 mt = event.moment_tensor.m_up_south_east()
175 ev_lat, ev_lon = event.effective_latlon
177 xx = num.trace(mt) / 3.
178 mc = num.matrix([[xx, 0., 0.], [0., xx, 0.], [0., 0., xx]])
179 mc = mt - mc
180 mc = mc / event.moment_tensor.scalar_moment() * \
181 moment_tensor.magnitude_to_moment(5.0)
182 m6 = tuple(moment_tensor.to6(mc))
183 symbol_size = 20.
184 m.gmt.psmeca(
185 S='%s%g' % ('d', symbol_size / gmtpy.cm),
186 in_rows=[(ev_lon, ev_lat, 10) + m6 + (1, 0, 0)],
187 M=True,
188 *m.jxyr)
190 elif isinstance(problem, RectangularProblem):
191 m.gmt.psxy(
192 in_rows=source.outline(cs='lonlat'),
193 L='+p2p,black',
194 W='1p,black',
195 G='black',
196 t=60,
197 *m.jxyr)
199 elif isinstance(problem, VLVDProblem):
200 ev_lat, ev_lon = event.effective_latlon
201 dV = abs(source.volume_change)
202 sphere_radius = num.cbrt(dV / (4./3.*num.pi))
204 volcanic_circle = [
205 ev_lon,
206 ev_lat,
207 '%fe' % sphere_radius
208 ]
209 m.gmt.psxy(
210 S='E-',
211 in_rows=[volcanic_circle],
212 W='1p,black',
213 G='orange3',
214 *m.jxyr)
216 return (item, m)
218 ifig = 0
219 for vertical in (False, True):
220 for gnss_target, result in zip(problem.gnss_targets, results):
221 yield plot_gnss(gnss_target, result, ifig, vertical)
222 ifig += 1
225class GNSSStationDistribution(StationDistributionPlot):
226 ''' Polar plot showing GNSS station distribution and weight '''
227 name = 'gnss_station_distribution'
229 def make(self, environ):
230 cm = environ.get_plot_collection_manager()
231 mpl_init(fontsize=self.font_size)
233 history = environ.get_history(subset='harvest')
234 problem = environ.get_problem()
235 dataset = environ.get_dataset()
237 cm.create_group_mpl(
238 self,
239 self.draw_figures(problem, dataset, history),
240 title=u'GNSS Station Distribution',
241 section='checks',
242 feather_icon='target',
243 description=u'''
244Plots showing the GNSS station distribution and their weight.
246This polar plot visualises the station distribution in distance and azimuth,
247the marker's size is scaled to the stations weight (mean of spatial
248components).
249''')
251 def draw_figures(self, problem, dataset, history):
253 event = problem.base_source
254 targets = problem.gnss_targets
256 for target in targets:
257 target.set_dataset(dataset)
258 comp_weights = target.component_weights()[0]
260 ws_n = comp_weights[:, 0::3] / comp_weights.max()
261 ws_e = comp_weights[:, 1::3] / comp_weights.max()
262 ws_u = comp_weights[:, 2::3] / comp_weights.max()
263 ws_e = num.array(ws_e[0]).flatten()
264 ws_n = num.array(ws_n[0]).flatten()
265 ws_u = num.array(ws_u[0]).flatten()
267 if ws_n.size == 0:
268 continue
270 distances = target.distance_to(event)
271 azimuths = od.azibazi_numpy(
272 num.array(event.effective_lat)[num.newaxis],
273 num.array(event.effective_lon)[num.newaxis],
274 target.get_latlon()[:, 0],
275 target.get_latlon()[:, 1])[0]
276 labels = target.station_names
278 item = PlotItem(name='station_distribution-N-%s' % target.path)
279 fig, ax, legend = self.plot_station_distribution(
280 azimuths, distances, ws_n, labels,
281 legend_title='Weight, N components')
283 yield (item, fig)
285 item = PlotItem(name='station_distribution-E-%s' % target.path)
286 fig, ax, legend = self.plot_station_distribution(
287 azimuths, distances, ws_e, labels,
288 legend_title='Weight, E components')
290 yield (item, fig)
292 item = PlotItem(name='station_distribution-U-%s' % target.path)
293 fig, ax, legend = self.plot_station_distribution(
294 azimuths, distances, ws_u, labels,
295 legend_title='Weight, U components')
297 yield (item, fig)
300def get_plot_classes():
301 return [
302 GNSSTargetMisfitPlot,
303 GNSSStationDistribution
304 ]