Merge pull request #9980 from gem/plot_rupture

ARISTOTLE: plot assets, rupture and hypocenter together...

openquake/calculators/event_based.py

@@ -315,8 +315,8 @@ def starmap_from_rups(func, oq, full_lt, sitecol, dstore, save_tmp=None):
     set_mags(oq, dstore)
     rups = dstore['ruptures'][:]
     logging.info('Reading {:_d} ruptures'.format(len(rups)))
-    logging.info('Affected sites ~%.0f per rupture, max=%.0f', rups['nsites'].mean(),
-                 rups['nsites'].max())
+    logging.info('Affected sites ~%.0f per rupture, max=%.0f',
+                 rups['nsites'].mean(), rups['nsites'].max())
     allproxies = [RuptureProxy(rec) for rec in rups]
     if "station_data" in oq.inputs:
         trt = full_lt.trts[0]
@@ -764,8 +764,9 @@ def save_avg_gmf(self):
         # make avg_gmf plots only if running via the webui
         if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE':
             imts = list(self.oqparam.imtls)
+            ex = Extractor(self.datastore.calc_id)
             for imt in imts:
-                plt = plot_avg_gmf(Extractor(self.datastore.calc_id), imt)
+                plt = plot_avg_gmf(ex, imt)
                 bio = io.BytesIO()
                 plt.savefig(bio, format='png', bbox_inches='tight')
                 fig_path = f'png/avg_gmf-{imt}.png'

openquake/calculators/postproc/plots.py

@@ -21,6 +21,7 @@
 from shapely.geometry import MultiPolygon
 from openquake.commonlib import readinput, datastore
 from openquake.hmtk.plotting.patch import PolygonPatch
+from openquake.calculators.getters import get_ebrupture
 
 
 def import_plt():
@@ -60,8 +61,8 @@ def get_country_iso_codes(calc_id, assetcol):
 
 def plot_avg_gmf(ex, imt):
     plt = import_plt()
-    fig = plt.figure()
-    ax = fig.add_subplot(1, 1, 1)
+    _fig, ax = plt.subplots(figsize=(10, 10))
+    ax.set_aspect('equal')
     ax.grid(True)
     ax.set_xlabel('Lon')
     ax.set_ylabel('Lat')
@@ -93,6 +94,34 @@ def plot_avg_gmf(ex, imt):
     return plt
 
 
+def add_rupture(ax, dstore, rup_id=0):
+    ebr = get_ebrupture(dstore, rup_id)
+    rup = ebr.rupture
+    poly = rup.surface.mesh.get_convex_hull()
+    min_x, min_y, max_x, max_y = poly.get_bbox()
+    ax.fill(poly.lons, poly.lats, alpha=.5, color='purple',
+            label='Rupture')
+    ax.plot(rup.hypocenter.x, rup.hypocenter.y, marker='*',
+            color='orange', label='Hypocenter', alpha=.5,
+            linestyle='', markersize=8)
+    return ax, min_x, min_y, max_x, max_y
+
+
+def plot_rupture(dstore):
+    plt = import_plt()
+    _fig, ax = plt.subplots(figsize=(10, 10))
+    ax.set_aspect('equal')
+    ax.grid(True)
+    # assuming there is only 1 rupture, so rup_id=0
+    ax, min_x, min_y, max_x, max_y = add_rupture(ax, dstore, rup_id=0)
+    ax = add_borders(ax)
+    BUF_ANGLE = 4
+    ax.set_xlim(min_x - BUF_ANGLE, max_x + BUF_ANGLE)
+    ax.set_ylim(min_y - BUF_ANGLE, max_y + BUF_ANGLE)
+    ax.legend()
+    return plt
+
+
 def get_assetcol(calc_id):
     try:
         dstore = datastore.read(calc_id)

openquake/commands/plot.py

@@ -32,7 +32,7 @@
 from openquake.calculators.extract import (
     Extractor, WebExtractor, clusterize)
 from openquake.calculators.postproc.plots import (
-    plot_avg_gmf, import_plt, add_borders)
+    plot_avg_gmf, import_plt, add_borders, plot_rupture)
 from openquake.calculators.postproc.aelo_plots import (
     plot_mean_hcurves_rtgm, plot_disagg_by_src, plot_governing_mce)
 
@@ -917,42 +917,12 @@ def plot_point_sources(srcs, ax, min_x, max_x, min_y, max_y):
     return min_x, max_x, min_y, max_y
 
 
-def make_figure_sources(extractors, what):
-    """
-    $ oq plot "sources?source_id=xxx"
-    $ oq plot "sources?code=N&code=F"
-    $ oq plot "sources?exclude=A"
-    """
-    # NB: matplotlib is imported inside since it is a costly import
-    plt = import_plt()
-    [ex] = extractors
-    dstore = ex.dstore
-    kwargs = what.split('?')[1]
-    if kwargs and 'exclude' in kwargs:
-        excluded_codes = [code.encode('utf8')
-                          for code in parse_qs(kwargs)['exclude']]
-    else:
-        excluded_codes = []
-    if kwargs and 'source_id' in kwargs:
-        src_ids = list(parse_qs(kwargs)['source_id'])
-    else:
-        src_ids = []
-    if kwargs and 'code' in kwargs:
-        codes = [code.encode('utf8') for code in parse_qs(kwargs)['code']
-                 if code.encode('utf8') not in excluded_codes]
-    else:
-        codes = []
+def plot_sources(srcs, ax, min_x, max_x, min_y, max_y):
     PLOTTABLE_CODES = (
         b'N', b'M', b'X', b'S', b'C', b'P', b'p', b'P', b'F', b'A', b'K')
-    print('Reading sources...')
-    csm = dstore['_csm']
-    srcs = filter_sources(csm, src_ids, codes, excluded_codes)
-    assert srcs, ('All sources were filtered out')
-    _fig, ax = plt.subplots()
     ax.set_aspect('equal')
     ax.grid(True)
     print(f'Plotting {len(srcs)} sources...')
-    min_x, max_x, min_y, max_y = (180, -180, 90, -90)
     any_sources_were_plotted = False
     # NonParametricSource
     np_sources = [src for src in srcs if src.code == b'N']
@@ -1006,6 +976,7 @@ def make_figure_sources(extractors, what):
     # MultiFaultSource
     mf_sources = [src for src in srcs if src.code == b'F']
     if mf_sources:
+        src_ids = [src.source_id for src in srcs]
         min_x, max_x, min_y, max_y = plot_multi_fault_sources(
             mf_sources, src_ids, ax, min_x, max_x, min_y, max_y)
         any_sources_were_plotted = True
@@ -1015,6 +986,42 @@ def make_figure_sources(extractors, what):
         print(f'Plotting the following sources is not'
               f'implemented yet: {unplottable}')
     assert any_sources_were_plotted, 'No sources were plotted'
+    return min_x, max_x, min_y, max_y
+
+
+def make_figure_sources(extractors, what):
+    """
+    $ oq plot "sources?source_id=xxx"
+    $ oq plot "sources?code=N&code=F"
+    $ oq plot "sources?exclude=A"
+    """
+    # NB: matplotlib is imported inside since it is a costly import
+    plt = import_plt()
+    [ex] = extractors
+    dstore = ex.dstore
+    kwargs = what.split('?')[1]
+    if kwargs and 'exclude' in kwargs:
+        excluded_codes = [code.encode('utf8')
+                          for code in parse_qs(kwargs)['exclude']]
+    else:
+        excluded_codes = []
+    if kwargs and 'source_id' in kwargs:
+        src_ids = list(parse_qs(kwargs)['source_id'])
+    else:
+        src_ids = []
+    if kwargs and 'code' in kwargs:
+        codes = [code.encode('utf8') for code in parse_qs(kwargs)['code']
+                 if code.encode('utf8') not in excluded_codes]
+    else:
+        codes = []
+    print('Reading sources...')
+    csm = dstore['_csm']
+    srcs = filter_sources(csm, src_ids, codes, excluded_codes)
+    assert srcs, ('All sources were filtered out')
+    _fig, ax = plt.subplots()
+    min_x, max_x, min_y, max_y = (180, -180, 90, -90)
+    min_x, max_x, min_y, max_y = plot_sources(
+        srcs, ax, min_x, max_x, min_y, max_y)
     print('Plotting mosaic borders...')
     ax = add_borders(ax, readinput.read_mosaic_df, buffer=0.)
     ax.set_xlim(min_x - ZOOM_MARGIN, max_x + ZOOM_MARGIN)
@@ -1026,6 +1033,16 @@ def make_figure_sources(extractors, what):
     return plt
 
 
+def make_figure_rupture(extractors, what):
+    """
+    $ oq plot "rupture?"
+    """
+    # NB: matplotlib is imported inside since it is a costly import
+    [ex] = extractors
+    dstore = ex.dstore
+    return plot_rupture(dstore)
+
+
 def plot_wkt(wkt_string):
     """
     Plot a WKT string describing a polygon

openquake/commands/plot_assets.py

@@ -16,13 +16,15 @@
 # You should have received a copy of the GNU Affero General Public License
 # along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.
 
+import os
 import numpy
 import shapely
 import logging
 from openquake.commonlib import datastore
 from openquake.hazardlib.geo.utils import cross_idl, get_bbox
 from openquake.calculators.postproc.plots import (
     add_borders, get_assetcol, get_country_iso_codes)
+from openquake.calculators.postproc.plots import add_rupture
 
 
 def main(calc_id: int = -1, site_model=False,
@@ -43,12 +45,12 @@ def main(calc_id: int = -1, site_model=False,
         region = None
     sitecol = dstore['sitecol']
     assetcol = get_assetcol(calc_id)
-    fig = p.figure()
-    ax = fig.add_subplot(111)
+    _fig, ax = p.subplots(figsize=(10, 10))
     if region:
         region_geom = shapely.wkt.loads(region)
         pp = PolygonPatch(region_geom, alpha=0.1)
         ax.add_patch(pp)
+    ax.set_aspect('equal')
     ax.grid(True)
     if assets_only:
         markersize_site_model = markersize_assets = 5
@@ -74,14 +76,18 @@ def main(calc_id: int = -1, site_model=False,
             p.scatter(disc['lon'], disc['lat'], marker='x', color='red',
                       label='discarded', s=markersize_discarded)
     if oq.rupture_xml or oq.rupture_dict:
-        # TODO: we may want to plot the rupture boundaries instead of just the
-        # hypocenter
         rec = dstore['ruptures'][0]
         lon, lat, _dep = rec['hypo']
         xlon, xlat = [lon], [lat]
         dist = sitecol.get_cdist(rec)
         print('rupture(%s, %s), dist=%s' % (lon, lat, dist))
-        p.scatter(xlon, xlat, marker='o', color='red', label='hypocenter')
+        if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE':
+            # assuming there is only 1 rupture, so rup_id=0
+            ax, _min_x, _min_y, _max_x, _max_y = add_rupture(
+                ax, dstore, rup_id=0)
+        else:
+            p.scatter(xlon, xlat, marker='*', color='orange',
+                      label='hypocenter', alpha=.5)
     else:
         xlon, xlat = [], []
 
@@ -92,9 +98,9 @@ def main(calc_id: int = -1, site_model=False,
     else:
         minx, miny, maxx, maxy = get_bbox(
             assetcol['lon'], assetcol['lat'], xlon, xlat)
-    w, h = maxx - minx, maxy - miny
-    ax.set_xlim(minx - 0.2 * w, maxx + 0.2 * w)
-    ax.set_ylim(miny - 0.2 * h, maxy + 0.2 * h)
+    BUF_ANGLE = 1
+    ax.set_xlim(minx - BUF_ANGLE, maxx + BUF_ANGLE)
+    ax.set_ylim(miny - BUF_ANGLE, maxy + BUF_ANGLE)
 
     country_iso_codes = get_country_iso_codes(calc_id, assetcol)
     if country_iso_codes is not None:

openquake/server/templates/engine/get_outputs_aristotle.html

@@ -33,7 +33,7 @@ <h2>Outputs from calculation {{ calc_id }}: {{ description }}{% if local_timesta
               {% if assets %}
               <div class="my-pngs">
                 <a href="{{ oq_engine_server_url }}/v1/calc/{{ calc_id }}/download_png/assets.png" class="btn btn-sm">
-                  Show assets</a>
+                  Show assets and rupture</a>
               </div>
               {% endif %}
             </div>

openquake/server/views.py

@@ -1566,7 +1566,8 @@ def web_engine_get_outputs_aristotle(request, calc_id):
                        time_job_after_event=time_job_after_event_str,
                        size_mb=size_mb, losses=losses,
                        losses_header=losses_header,
-                       avg_gmf=avg_gmf, assets=assets, warnings=warnings))
+                       avg_gmf=avg_gmf, assets=assets,
+                       warnings=warnings))
 
 
 @cross_domain_ajax