A probabilistic earthquake source inversion framework. Designed and crafted in Mordor. https://pyrocko.org/grond/
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

408 lines
15 KiB

  1. import logging
  2. import numpy as num
  3. from matplotlib import cm, gridspec
  4. from grond.plot.config import PlotConfig
  5. from grond.plot.collection import PlotItem
  6. from matplotlib import pyplot as plt
  7. from matplotlib.ticker import MaxNLocator
  8. from matplotlib import patches
  9. from pyrocko.guts import Tuple, Float, String, Int, Bool, StringChoice
  10. logger = logging.getLogger('grond.targets.satellite.plot')
  11. km = 1e3
  12. d2r = num.pi/180.
  13. guts_prefix = 'grond'
  14. def scale_axes(axis, scale, offset=0., suffix=''):
  15. from matplotlib.ticker import ScalarFormatter
  16. class FormatScaled(ScalarFormatter):
  17. @staticmethod
  18. def __call__(value, pos):
  19. return '{:,.1f}{:}'.format((offset + value) * scale, suffix)\
  20. .replace(',', ' ')
  21. axis.set_major_formatter(FormatScaled())
  22. class SatelliteTargetDisplacement(PlotConfig):
  23. ''' Maps showing surface displacements from satellite and modelled data '''
  24. name = 'satellite'
  25. dpi = Int.T(
  26. default=250)
  27. size_cm = Tuple.T(
  28. 2, Float.T(),
  29. default=(22., 12.))
  30. colormap = String.T(
  31. default='RdBu',
  32. help='Colormap for the surface displacements')
  33. relative_coordinates = Bool.T(
  34. default=False,
  35. help='Show relative coordinates, initial location centered at 0N, 0E')
  36. fit = StringChoice.T(
  37. default='best', choices=['best', 'mean'],
  38. help='Show the \'best\' or \'mean\' fits and source model from the'
  39. ' ensamble.')
  40. def make(self, environ):
  41. cm = environ.get_plot_collection_manager()
  42. history = environ.get_history(subset='harvest')
  43. optimiser = environ.get_optimiser()
  44. ds = environ.get_dataset()
  45. environ.setup_modelling()
  46. cm.create_group_mpl(
  47. self,
  48. self.draw_static_fits(ds, history, optimiser),
  49. title=u'InSAR Displacements',
  50. section='fits',
  51. feather_icon='navigation',
  52. description=u'''
  53. Maps showing subsampled surface displacements as observed, modelled and the
  54. residual (observed minus modelled).
  55. The displacement values predicted by the orbit-ambiguity ramps are added to the
  56. modelled displacements (middle panels). The color shows the LOS displacement
  57. values associated with, and the extent of, every quadtree box. The light grey
  58. dots show the focal point of pixels combined in the quadtree box. This point
  59. corresponds to the position of the modelled data point.
  60. The large dark grey dot shows the reference source position. The grey filled
  61. box shows the surface projection of the modelled source, with the thick-lined
  62. edge marking the upper fault edge. Complete data extent is shown.
  63. ''')
  64. def draw_static_fits(self, ds, history, optimiser, closeup=False):
  65. from pyrocko.orthodrome import latlon_to_ne_numpy
  66. problem = history.problem
  67. sat_targets = problem.satellite_targets
  68. for target in sat_targets:
  69. target.set_dataset(ds)
  70. if self.fit == 'best':
  71. source = history.get_best_source()
  72. model = history.get_best_model()
  73. elif self.fit == 'mean':
  74. source = history.get_mean_source()
  75. model = history.get_mean_model()
  76. results = problem.evaluate(model, targets=sat_targets)
  77. def initAxes(ax, scene, title, last_axes=False):
  78. ax.set_title(title)
  79. ax.tick_params(length=2)
  80. if scene.frame.isMeter():
  81. ax.set_xlabel('Easting [km]')
  82. scale_x = dict(scale=1./km)
  83. scale_y = dict(scale=1./km)
  84. if not self.relative_coordinates:
  85. import utm
  86. utm_E, utm_N, utm_zone, utm_zone_letter =\
  87. utm.from_latlon(source.effective_lat,
  88. source.effective_lon)
  89. scale_x['offset'] = utm_E
  90. scale_y['offset'] = utm_N
  91. if last_axes:
  92. ax.text(0.975, 0.025,
  93. 'UTM Zone %d%s' % (utm_zone, utm_zone_letter),
  94. va='bottom', ha='right',
  95. fontsize=8, alpha=.7,
  96. transform=ax.transAxes)
  97. ax.set_aspect('equal')
  98. elif scene.frame.isDegree():
  99. scale_x = dict(scale=1., suffix='°')
  100. scale_y = dict(scale=1., suffix='°')
  101. if not self.relative_coordinates:
  102. scale_x['offset'] = source.effective_lon
  103. scale_y['offset'] = source.effective_lat
  104. ax.set_aspect(1./num.cos(source.effective_lat*d2r))
  105. nticks_lon = 4 if abs(scene.frame.llLon) >= 100 else 5
  106. ax.xaxis.set_major_locator(MaxNLocator(nticks_lon))
  107. ax.yaxis.set_major_locator(MaxNLocator(5))
  108. scale_axes(ax.get_xaxis(), **scale_x)
  109. scale_axes(ax.get_yaxis(), **scale_y)
  110. def drawSource(ax, scene):
  111. if scene.frame.isMeter():
  112. fn, fe = source.outline(cs='xy').T
  113. fn -= fn.mean()
  114. fe -= fe.mean()
  115. elif scene.frame.isDegree():
  116. fn, fe = source.outline(cs='latlon').T
  117. fn -= source.effective_lat
  118. fe -= source.effective_lon
  119. # source is centered
  120. ax.scatter(0., 0., color='black', s=3, alpha=.5, marker='o')
  121. ax.fill(fe, fn,
  122. edgecolor=(0., 0., 0.),
  123. facecolor=(.5, .5, .5), alpha=0.7)
  124. ax.plot(fe[0:2], fn[0:2], 'k', linewidth=1.3)
  125. def mapDisplacementGrid(displacements, scene):
  126. arr = num.full_like(scene.displacement, fill_value=num.nan)
  127. qt = scene.quadtree
  128. for syn_v, l in zip(displacements, qt.leaves):
  129. arr[l._slice_rows, l._slice_cols] = syn_v
  130. arr[scene.displacement_mask] = num.nan
  131. return arr
  132. def drawLeaves(ax, scene, offset_e=0., offset_n=0.):
  133. rects = scene.quadtree.getMPLRectangles()
  134. for r in rects:
  135. r.set_edgecolor((.4, .4, .4))
  136. r.set_linewidth(.5)
  137. r.set_facecolor('none')
  138. r.set_x(r.get_x() - offset_e)
  139. r.set_y(r.get_y() - offset_n)
  140. map(ax.add_artist, rects)
  141. ax.scatter(scene.quadtree.leaf_coordinates[:, 0] - offset_e,
  142. scene.quadtree.leaf_coordinates[:, 1] - offset_n,
  143. s=.25, c='black', alpha=.1)
  144. def addArrow(ax, scene):
  145. phi = num.nanmean(scene.phi)
  146. los_dx = num.cos(phi + num.pi) * .0625
  147. los_dy = num.sin(phi + num.pi) * .0625
  148. az_dx = num.cos(phi - num.pi/2) * .125
  149. az_dy = num.sin(phi - num.pi/2) * .125
  150. anchor_x = .9 if los_dx < 0 else .1
  151. anchor_y = .85 if los_dx < 0 else .975
  152. az_arrow = patches.FancyArrow(
  153. x=anchor_x-az_dx, y=anchor_y-az_dy,
  154. dx=az_dx, dy=az_dy,
  155. head_width=.025,
  156. alpha=.5, fc='k',
  157. head_starts_at_zero=False,
  158. length_includes_head=True,
  159. transform=ax.transAxes)
  160. los_arrow = patches.FancyArrow(
  161. x=anchor_x-az_dx/2, y=anchor_y-az_dy/2,
  162. dx=los_dx, dy=los_dy,
  163. head_width=.02,
  164. alpha=.5, fc='k',
  165. head_starts_at_zero=False,
  166. length_includes_head=True,
  167. transform=ax.transAxes)
  168. ax.add_artist(az_arrow)
  169. ax.add_artist(los_arrow)
  170. urE, urN, llE, llN = (0., 0., 0., 0.)
  171. for target in sat_targets:
  172. if target.scene.frame.isMeter():
  173. off_n, off_e = map(float, latlon_to_ne_numpy(
  174. target.scene.frame.llLat, target.scene.frame.llLon,
  175. source.effective_lat, source.effective_lon))
  176. if target.scene.frame.isDegree():
  177. off_n = source.effective_lat - target.scene.frame.llLat
  178. off_e = source.effective_lon - target.scene.frame.llLon
  179. turE, turN, tllE, tllN = zip(
  180. *[(l.gridE.max()-off_e,
  181. l.gridN.max()-off_n,
  182. l.gridE.min()-off_e,
  183. l.gridN.min()-off_n)
  184. for l in target.scene.quadtree.leaves])
  185. turE, turN = map(max, (turE, turN))
  186. tllE, tllN = map(min, (tllE, tllN))
  187. urE, urN = map(max, ((turE, urE), (urN, turN)))
  188. llE, llN = map(min, ((tllE, llE), (llN, tllN)))
  189. def generate_plot(sat_target, result, ifig):
  190. scene = sat_target.scene
  191. fig = plt.figure()
  192. fig.set_size_inches(*self.size_inch)
  193. gs = gridspec.GridSpec(
  194. 2, 3,
  195. wspace=.15, hspace=.2,
  196. left=.1, right=.975, top=.95,
  197. height_ratios=[12, 1])
  198. item = PlotItem(
  199. name='fig_%i' % ifig,
  200. attributes={'targets': [sat_target.path]},
  201. title=u'Satellite Surface Displacements - %s'
  202. % scene.meta.scene_title,
  203. description=u'''
  204. Surface displacements derived from satellite data.
  205. (Left) the input data, (center) the modelled
  206. data and (right) the model residual.
  207. '''.format(meta=scene.meta))
  208. stat_obs = result.statics_obs
  209. stat_syn = result.statics_syn['displacement.los']
  210. res = stat_obs - stat_syn
  211. if scene.frame.isMeter():
  212. offset_n, offset_e = map(float, latlon_to_ne_numpy(
  213. scene.frame.llLat, scene.frame.llLon,
  214. source.effective_lat, source.effective_lon))
  215. elif scene.frame.isDegree():
  216. offset_n = source.effective_lat - scene.frame.llLat
  217. offset_e = source.effective_lon - scene.frame.llLon
  218. im_extent = (scene.frame.E.min() - offset_e,
  219. scene.frame.E.max() - offset_e,
  220. scene.frame.N.min() - offset_n,
  221. scene.frame.N.max() - offset_n)
  222. abs_displ = num.abs([stat_obs.min(), stat_obs.max(),
  223. stat_syn.min(), stat_syn.max(),
  224. res.min(), res.max()]).max()
  225. cmw = cm.ScalarMappable(cmap=self.colormap)
  226. cmw.set_clim(vmin=-abs_displ, vmax=abs_displ)
  227. cmw.set_array(stat_obs)
  228. axes = [fig.add_subplot(gs[0, 0]),
  229. fig.add_subplot(gs[0, 1]),
  230. fig.add_subplot(gs[0, 2])]
  231. ax = axes[0]
  232. ax.imshow(mapDisplacementGrid(stat_obs, scene),
  233. extent=im_extent, cmap=self.colormap,
  234. vmin=-abs_displ, vmax=abs_displ,
  235. origin='lower')
  236. drawLeaves(ax, scene, offset_e, offset_n)
  237. drawSource(ax, scene)
  238. addArrow(ax, scene)
  239. initAxes(ax, scene, 'Observed')
  240. ax.text(.025, .025, 'Scene ID: %s' % scene.meta.scene_id,
  241. fontsize=8, alpha=.7,
  242. va='bottom', transform=ax.transAxes)
  243. if scene.frame.isMeter():
  244. ax.set_ylabel('Northing [km]')
  245. ax = axes[1]
  246. ax.imshow(mapDisplacementGrid(stat_syn, scene),
  247. extent=im_extent, cmap=self.colormap,
  248. vmin=-abs_displ, vmax=abs_displ,
  249. origin='lower')
  250. drawLeaves(ax, scene, offset_e, offset_n)
  251. drawSource(ax, scene)
  252. addArrow(ax, scene)
  253. initAxes(ax, scene, 'Model')
  254. ax.get_yaxis().set_visible(False)
  255. ax = axes[2]
  256. ax.imshow(mapDisplacementGrid(res, scene),
  257. extent=im_extent, cmap=self.colormap,
  258. vmin=-abs_displ, vmax=abs_displ,
  259. origin='lower')
  260. drawLeaves(ax, scene, offset_e, offset_n)
  261. drawSource(ax, scene)
  262. addArrow(ax, scene)
  263. initAxes(ax, scene, 'Residual', last_axes=True)
  264. ax.get_yaxis().set_visible(False)
  265. for ax in axes:
  266. ax.set_xlim(llE, urE)
  267. ax.set_ylim(llN, urN)
  268. if closeup:
  269. if scene.frame.isMeter():
  270. fn, fe = source.outline(cs='xy').T
  271. elif scene.frame.isDegree():
  272. fn, fe = source.outline(cs='latlon').T
  273. fn -= source.effective_lat
  274. fe -= source.effective_lon
  275. if fn.size > 1:
  276. off_n = (fn[0] + fn[1]) / 2
  277. off_e = (fe[0] + fe[1]) / 2
  278. else:
  279. off_n = fn[0]
  280. off_e = fe[0]
  281. fault_size = 2*num.sqrt(max(abs(fn-off_n))**2
  282. + max(abs(fe-off_e))**2)
  283. fault_size *= self.map_scale
  284. if fault_size == 0.0:
  285. extent = (scene.frame.N[-1] + scene.frame.E[-1]) / 2
  286. fault_size = extent * .25
  287. for ax in axes:
  288. ax.set_xlim(-fault_size/2 + off_e, fault_size/2 + off_e)
  289. ax.set_ylim(-fault_size/2 + off_n, fault_size/2 + off_n)
  290. cax = fig.add_subplot(gs[1, :])
  291. cbar = fig.colorbar(cmw, cax=cax, orientation='horizontal',
  292. aspect=20, use_gridspec=True)
  293. cbar.set_label('LOS Displacement [m]')
  294. return (item, fig)
  295. for ifig, (sat_target, result) in enumerate(zip(sat_targets, results)):
  296. yield generate_plot(sat_target, result, ifig)
  297. class SatelliteTargetDisplacementCloseup(SatelliteTargetDisplacement):
  298. ''' Close-up of satellite surface displacements and modelled data. '''
  299. name = 'satellite_closeup'
  300. map_scale = Float.T(
  301. default=2.,
  302. help='Scale the map surroundings, larger value zooms out.')
  303. def make(self, environ):
  304. cm = environ.get_plot_collection_manager()
  305. history = environ.get_history(subset='harvest')
  306. optimiser = environ.get_optimiser()
  307. ds = environ.get_dataset()
  308. environ.setup_modelling()
  309. cm.create_group_mpl(
  310. self,
  311. self.draw_static_fits(ds, history, optimiser, closeup=True),
  312. title=u'InSAR Displacements (Closeup)',
  313. section='fits',
  314. feather_icon='zoom-in',
  315. description=u'''
  316. Maps showing subsampled surface displacements as observed, modelled and the
  317. residual (observed minus modelled).
  318. The displacement values predicted by the orbit-ambiguity ramps are added to the
  319. modelled displacements (middle panels). The color shows the LOS displacement
  320. values associated with, and the extent of, every quadtree box. The light grey
  321. dots show the focal point of pixels combined in the quadtree box. This point
  322. corresponds to the position of the modelled data point.
  323. The large dark grey dot shows the reference source position. The grey filled
  324. box shows the surface projection of the modelled source, with the thick-lined
  325. edge marking the upper fault edge. Map is focused around the fault's extent.
  326. ''')
  327. def get_plot_classes():
  328. return [SatelliteTargetDisplacement, SatelliteTargetDisplacementCloseup]