Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/gnss_campaign/plot.py: 22%

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import logging 

5import numpy as num 

6from pyrocko.model import gnss 

7 

8from pyrocko.plot import automap, mpl_init 

9from pyrocko import orthodrome as od 

10 

11from grond.plot.config import PlotConfig 

12from grond.plot.collection import PlotItem 

13from grond.problems import CMTProblem, RectangularProblem, \ 

14 VLVDProblem 

15 

16from ..plot import StationDistributionPlot 

17 

18import copy 

19from pyrocko.guts import Tuple, Float, Bool 

20 

21guts_prefix = 'grond' 

22km = 1e3 

23 

24logger = logging.getLogger('grond.targets.gnss_campaign.plot') 

25 

26 

27class GNSSTargetMisfitPlot(PlotConfig): 

28 ''' Maps showing horizontal surface displacements 

29 of a GNSS campaign and model ''' 

30 

31 name = 'gnss' 

32 

33 size_cm = Tuple.T( 

34 2, Float.T(), 

35 default=(30., 30.), 

36 help='width and length of the figure in cm') 

37 show_topo = Bool.T( 

38 default=False, 

39 help='show topography') 

40 show_grid = Bool.T( 

41 default=True, 

42 help='show the lat/lon grid') 

43 show_rivers = Bool.T( 

44 default=True, 

45 help='show rivers on the map') 

46 radius = Float.T( 

47 optional=True, 

48 help='radius of the map around campaign center lat/lon') 

49 

50 def make(self, environ): 

51 cm = environ.get_plot_collection_manager() 

52 history = environ.get_history(subset='harvest') 

53 optimiser = environ.get_optimiser() 

54 ds = environ.get_dataset() 

55 

56 environ.setup_modelling() 

57 

58 cm.create_group_automap( 

59 self, 

60 self.draw_gnss_fits(ds, history, optimiser), 

61 title=u'GNSS Displacements', 

62 section='fits', 

63 feather_icon='map', 

64 description=u''' 

65Maps showing station positions and statiom names of the GNSS targets. 

66 

67Arrows the observed surface displacements (black arrows) and synthetic 

68displacements (red arrows). The top plot shows the horizontal displacements and 

69the bottom plot the vertical displacements. The grey filled box shows the 

70surface projection of the modelled source, with the thick-lined edge marking 

71the upper fault edge. 

72''') 

73 

74 def draw_gnss_fits(self, ds, history, optimiser, vertical=False): 

75 problem = history.problem 

76 

77 gnss_targets = problem.gnss_targets 

78 for target in gnss_targets: 

79 target.set_dataset(ds) 

80 

81 xbest = history.get_best_model() 

82 source = history.get_best_source() 

83 

84 results = problem.evaluate( 

85 xbest, result_mode='full', targets=gnss_targets) 

86 

87 def plot_gnss(gnss_target, result, ifig, vertical=False): 

88 campaign = gnss_target.campaign 

89 item = PlotItem( 

90 name='fig_%i' % ifig, 

91 attributes={ 

92 'targets': gnss_target.path 

93 }, 

94 title=u'Static GNSS Surface Displacements - Campaign %s' 

95 % campaign.name, 

96 description=u''' 

97Static surface displacement from GNSS campaign %s (black vectors) and 

98displacements derived from best model (red). 

99''' % campaign.name) 

100 

101 event = source.pyrocko_event() 

102 locations = campaign.stations + [event] 

103 

104 lat, lon = od.geographic_midpoint_locations(locations) 

105 

106 if self.radius is None: 

107 coords = num.array([loc.effective_latlon for loc in locations]) 

108 radius = od.distance_accurate50m_numpy( 

109 lat[num.newaxis], lon[num.newaxis], 

110 coords[:, 0].max(), coords[:, 1]).max() 

111 radius *= 1.1 

112 

113 if radius < 30.*km: 

114 logger.warning( 

115 'Radius of GNSS campaign %s too small, defaulting' 

116 ' to 30 km' % campaign.name) 

117 radius = 30*km 

118 

119 model_camp = gnss.GNSSCampaign( 

120 stations=copy.deepcopy(campaign.stations), 

121 name='grond model') 

122 for ista, sta in enumerate(model_camp.stations): 

123 sta.north.shift = result.statics_syn['displacement.n'][ista] 

124 sta.north.sigma = 0. 

125 

126 sta.east.shift = result.statics_syn['displacement.e'][ista] 

127 sta.east.sigma = 0. 

128 

129 if sta.up: 

130 sta.up.shift = -result.statics_syn['displacement.d'][ista] 

131 sta.up.sigma = 0. 

132 

133 m = automap.Map( 

134 width=self.size_cm[0], 

135 height=self.size_cm[1], 

136 lat=lat, 

137 lon=lon, 

138 radius=radius, 

139 show_topo=self.show_topo, 

140 show_grid=self.show_grid, 

141 show_rivers=self.show_rivers, 

142 color_wet=(216, 242, 254), 

143 color_dry=(238, 236, 230)) 

144 

145 all_stations = campaign.stations + model_camp.stations 

146 offset_scale = num.zeros(len(all_stations)) 

147 

148 for ista, sta in enumerate(all_stations): 

149 for comp in sta.components.values(): 

150 offset_scale[ista] += comp.shift 

151 offset_scale = num.sqrt(offset_scale**2).max() 

152 

153 m.add_gnss_campaign( 

154 campaign, 

155 psxy_style={ 

156 'G': 'black', 

157 'W': '0.8p,black', 

158 }, 

159 offset_scale=offset_scale, 

160 vertical=vertical) 

161 

162 m.add_gnss_campaign( 

163 model_camp, 

164 psxy_style={ 

165 'G': 'red', 

166 'W': '0.8p,red', 

167 't': 30, 

168 }, 

169 offset_scale=offset_scale, 

170 vertical=vertical, 

171 labels=False) 

172 

173 if isinstance(problem, CMTProblem): 

174 from pyrocko import moment_tensor 

175 from pyrocko.plot import gmtpy 

176 

177 mt = event.moment_tensor.m_up_south_east() 

178 ev_lat, ev_lon = event.effective_latlon 

179 

180 xx = num.trace(mt) / 3. 

181 mc = num.matrix([[xx, 0., 0.], [0., xx, 0.], [0., 0., xx]]) 

182 mc = mt - mc 

183 mc = mc / event.moment_tensor.scalar_moment() * \ 

184 moment_tensor.magnitude_to_moment(5.0) 

185 m6 = tuple(moment_tensor.to6(mc)) 

186 symbol_size = 20. 

187 m.gmt.psmeca( 

188 S='%s%g' % ('d', symbol_size / gmtpy.cm), 

189 in_rows=[(ev_lon, ev_lat, 10) + m6 + (1, 0, 0)], 

190 M=True, 

191 *m.jxyr) 

192 

193 elif isinstance(problem, RectangularProblem): 

194 m.gmt.psxy( 

195 in_rows=source.outline(cs='lonlat'), 

196 L='+p2p,black', 

197 W='1p,black', 

198 G='black', 

199 t=60, 

200 *m.jxyr) 

201 

202 elif isinstance(problem, VLVDProblem): 

203 ev_lat, ev_lon = event.effective_latlon 

204 dV = abs(source.volume_change) 

205 sphere_radius = num.cbrt(dV / (4./3.*num.pi)) 

206 

207 volcanic_circle = [ 

208 ev_lon, 

209 ev_lat, 

210 '%fe' % sphere_radius 

211 ] 

212 m.gmt.psxy( 

213 S='E-', 

214 in_rows=[volcanic_circle], 

215 W='1p,black', 

216 G='orange3', 

217 *m.jxyr) 

218 

219 return (item, m) 

220 

221 ifig = 0 

222 for vertical in (False, True): 

223 for gnss_target, result in zip(problem.gnss_targets, results): 

224 yield plot_gnss(gnss_target, result, ifig, vertical) 

225 ifig += 1 

226 

227 

228class GNSSStationDistribution(StationDistributionPlot): 

229 ''' Polar plot showing GNSS station distribution and weight ''' 

230 name = 'gnss_station_distribution' 

231 

232 def make(self, environ): 

233 cm = environ.get_plot_collection_manager() 

234 mpl_init(fontsize=self.font_size) 

235 

236 history = environ.get_history(subset='harvest') 

237 problem = environ.get_problem() 

238 dataset = environ.get_dataset() 

239 

240 cm.create_group_mpl( 

241 self, 

242 self.draw_figures(problem, dataset, history), 

243 title=u'GNSS Station Distribution', 

244 section='checks', 

245 feather_icon='target', 

246 description=u''' 

247Plots showing the GNSS station distribution and their weight. 

248 

249This polar plot visualises the station distribution in distance and azimuth, 

250the marker's size is scaled to the stations weight (mean of spatial 

251components). 

252''') 

253 

254 def draw_figures(self, problem, dataset, history): 

255 

256 event = problem.base_source 

257 targets = problem.gnss_targets 

258 

259 for target in targets: 

260 target.set_dataset(dataset) 

261 comp_weights = target.component_weights()[0] 

262 

263 ws_n = comp_weights[:, 0::3] / comp_weights.max() 

264 ws_e = comp_weights[:, 1::3] / comp_weights.max() 

265 ws_u = comp_weights[:, 2::3] / comp_weights.max() 

266 ws_e = num.array(ws_e[0]).flatten() 

267 ws_n = num.array(ws_n[0]).flatten() 

268 ws_u = num.array(ws_u[0]).flatten() 

269 

270 if ws_n.size == 0: 

271 continue 

272 

273 distances = target.distance_to(event) 

274 azimuths = od.azibazi_numpy( 

275 num.array(event.effective_lat)[num.newaxis], 

276 num.array(event.effective_lon)[num.newaxis], 

277 target.get_latlon()[:, 0], 

278 target.get_latlon()[:, 1])[0] 

279 labels = target.station_names 

280 

281 item = PlotItem(name='station_distribution-N-%s' % target.path) 

282 fig, ax, legend = self.plot_station_distribution( 

283 azimuths, distances, ws_n, labels, 

284 legend_title='Weight, N components') 

285 

286 yield (item, fig) 

287 

288 item = PlotItem(name='station_distribution-E-%s' % target.path) 

289 fig, ax, legend = self.plot_station_distribution( 

290 azimuths, distances, ws_e, labels, 

291 legend_title='Weight, E components') 

292 

293 yield (item, fig) 

294 

295 item = PlotItem(name='station_distribution-U-%s' % target.path) 

296 fig, ax, legend = self.plot_station_distribution( 

297 azimuths, distances, ws_u, labels, 

298 legend_title='Weight, U components') 

299 

300 yield (item, fig) 

301 

302 

303def get_plot_classes(): 

304 return [ 

305 GNSSTargetMisfitPlot, 

306 GNSSStationDistribution 

307 ]